MATLAB課題 

D5 25 内藤 聡

 

% Program 1 --- Expression of transfer function ---

>> num=5*[1 3]; den=conv([1 6], [1 8]) ; printsys(num,den,'s');

分子を5(s+3) 分母を(s+6)(s+8)として伝達関数H(Z)を表示

num/den =

 

% Program 2 --- Transformation from frequency domain to space state domain ---

前問で使ったnumとdenを使って状態区間型から伝達関数型に変換する

>> [aa,bb,cc,dd]=tf2ss(num, den)

aa =

-14 -48

1 0

bb =

1

0

cc =

5 15

dd =

0

 

% Program 3 --- Expression of state space model ---

状態空間モデルで表現する

>> a=[-3 -4 -12; 1 0 0; 0 1 0]; b=[1 0 0]'; c=[1 3 2]; d=0; printsys(a,b,c,d);

a =

x1 x2 x3

x1 -3.00000 -4.00000 -12.00000

x2 1.00000 0 0

x3 0 1.00000 0

b =

u1

x1 1.00000

x2 0

x3 0

c =

x1 x2 x3

y1 1.00000 3.00000 2.00000

d =

u1

y1 0

 

% Program 4 --- Transformation from state space domain to frequency domain

状態空間型から伝達関数型へ変換

>> [numm, denn] = ss2tf(a,b,c,d,1);printsys(numm,denn,'s')

num/den =

 

% Program 5 --- Series and paralled connection of two models ---

4で求めた伝達関数をseriesで直列接続、parallelで並列接続

>> [nums,dens]=series(numm,denn,numm,denn); printsys(nums,dens,'s')

num/den =

 

>>[nump,denp]=parallel(numm,denn,numm,denn); printsys(nump,denp,'s')

num/den =

 

% Program 6 --- Closed loop transfer function for unit feedback system ---

4で求めた伝達関数の閉ループ伝達関数を求める

>> [numc,denc]=cloop(numm,denn);printsys(numc,denc,'s')

num/den =

 

% Program 7 --- Closed loop transfer function for unit feedback system ---

pro1の伝達関数とpro4の伝達関数の結合

>> [numf,denf] = feedback(num,den,numm,denn); printsys(numf, denf, 's')

num/den =

 

% Program 8 --- Reduction from subsystems, to single transfer function

各成分と行列q、入出力数を決定、伝達空間型から状態空間型へ変換

>>n1=1;d1=1;n2=1;d2=[11];n3=1;d3=[12];n4=1;d4=[13];n5=4;d5=1;n6=8;d6=1;n7=12;d7=1;nblocks=7;blkbuild

State model [a,b,c,d] of the block diagram has 7 inputs and 7 outputs

>> q=[2 -5 1 0 0; 3 2 -6 0 0; 4 3 -6 2 -7

5 3 0 0 0; 6 3 0 0 0; 7 4 0 0 0];

>> input=1; output=4; [aa,bb,cc,dd]=connect(a,b,c,d,q,input,output);

>> [num,den]=ss2tf(aa,bb,cc,dd); printsys(num,den,'s')

num/den =

 

% Program 9 --- Computing poles and zeros ---

状態空間型から伝達空間型へ変換。ゼロ点、極、利得型から伝達関数型へ

>> [z,p,k]=tf2zp(num,den)

z =

-3.0000

p =

-15.0000

-9.5311

-1.4689

k =

1.0000

>> [z,p,k]=ss2zp(aa,bb,cc,dd)

z =

-3.0000

p =

-15.0000

-9.5311

-1.4689

k =

1

 

% Program 10 --- Creating transfer function from poles and zeros ---

伝達関数型からゼロ点、極、利得型へ

>> [numcr,dencr]=zp2tf(z,p,k);printsys(numcr,dencr,'s');

num/den =

 

% Program 11 --- Step function and plot ---

伝達関数型からゼロ点、極、利得型へ変換

>> [num,den] = zp2tf([],-2+5*i -2+5*i,2); step(num,den)

 

% Program 12 --- Analysis of step response ---

伝達関数を定義して最終値定理を求める。また 立ち上がり時間と整定時間も求める

>> [num,den]=zp2tf([],[-1+3*i -1-3*i],3); finalv=polyval(num,0)/polyval(den,0)

finalv =

0.3000

>> [y,x,t]=step(num,den); [Y,k]=max(y);

>> peak=t(k)

peak =

1.0491

>> percent=100*(Y-finalv)/finalv

percent =

35.0914

% Compute rise time

>> n=1; while y(n)<0.1*finalv,n=n+1;end

>> m=1; while y(m)<0.9*finalv,m=m+1;end

>> rtime=t(m)-t(n)

rtime =

0.4417

% Compute settling time

>> l=length(t); while (y(l)>0.98*finalv)&(y(l)<1.02*finalv)

l=l-1; end

>> stime=t(l)

stime =

3.5337

>> pause

 

% Program 13 --- Meshplot of second-order system ---

グラフ化

>> clg

Warning: This function is obsolete and may be removed in future versions.

Use CLF instead.

> In /usr/matlab5/toolbox/matlab/graphics/clg.m at line 11

>> clf

>> t=[0:0.05:5]; number=10; y=zeros(length(t),number); n=1; while n<=number,

[num,den]=zp2tf([],[-n/4+2.5*i -n/4-2.5*i],(n/4)^2+8);

[y(l:length(t),n)x,tdumb]=step(num,den,t);

n=n+1;

end

clf

mesh(y')

title('Mesh Plot Showing Step Response for Twelve Pole Locations')

pause

%Result

??? In an assignment A(matrix,matrix) = B, the number of rows in B

and the number of elements in the A row index matrix must be the same.

 

% Program 14 --- Plot of time history reponse and phase plane ---

二次システムに対し行列を作成。グラフ化

clf

[num,den]=ord2(10,0.1); subplot(221)

step(num,den)

subplot(222)

impulse(num,den)

temp=fliplr(den); a=[0 1;-temp([1 1 0])]; b=[0;1]; c=[1 0]; d=0;

%error

??? Index into matrix is negative or zero. See release notes on changes to

logical indices.

 

[y,x,t]=step(a,b,c,d); subplot(223)

plot(x(:,1),x(:,2))

xlabel('x1-axis'),ylabel('x2-axis'),

title('x2 vs. x1 for step input'),

[y1,x1,t1]=impulse(a,b,c,d); subplot(224),

plot(x1(:,1),x1(:,2))

xlabel('x1-axis'),ylabel('x2-axis')

title('x2 vs. x1 for impulse input')

% Program 15 --- Root Locus plot ---

evansの根軌跡を求める

>> clf

>> num=1; den=conv(conv([1 0.5],[1 1]),[0.5 1]); [r,kvector]=rlocus(num,den);

>> v1=0.5; v2=2.5; h1=2.5; h2=0.2; axis([-h1 h2 -v1 v2]); plot(r,'-')

>> xlabel('Real Axis')

>> ylabel('Imaginary Axis')

>> hold on

>> plot(real(r(1,:)),imag(r(1,:)),'x')

>> plot([-h1 h2],[0 0],'w:')

>> plot([0 0],[-v1 v2],'w:')

>> hold off

>> dmp=0.707; wn=1:1:4; sgrid(dmp,wn)

% Program 16 --- Nyquist plot and bode plot---

ボード線図を求める

>> clf

>> [num,den]=zp2tf([],[-0.5 -1 -2],5); printsys(num,den,'s')

num/den =

 

>> [numgc,dengc]=cloop(num,den); printsys(numgc,dengc,'s')

num/den =

 

>> nyquist(num,den)

>> pause

>> clf

>> bode(num,den)

>> w=logspace(-1,1.100);

>> [mag,phase,w]=bode(num,den);[gainm,phasem,wc,wp]=margin(mag,phase,w)

gainm =

2.2623

 

phasem =

29.2403

 

wc =

1.8729

 

wp =

1.2369

>> pause; bode(numgc,dengc)

>> [maggc,phasegc,w]=bode(numgc,dengc); [mp,k]=max(maggc)

mp =

2.1785

 

k =

21

>> wp=w(k)

wp =

1.3693

>> n=1; while 20*log10(maggc(n))>=-3,n=n+1; end

>> bandwidth=w(n)

bandwidth =

1.9951

% Program 17 ---Nichols chart---

ニコラス線図を求める

num=1;

den=conv([1 10],conv([1 0.0001],[4 3]));

axis([-360 0 -40 80])

%[mag,phase,w]=nichols(num,den);

[mag1,phase1,w1]=nichols(100*num,den);

plot(phase,20*log10(mag),'-',phase1,20*log10(mag1),'--')

title('Nichols plot')

ngrid

% program18---contorol system design using phase compensator---

program

制御対象

位相補償器

閉ループ伝達関数

閉ループ伝達関数

制御前後のステップ応答

%program19---control system design using PID compensator---

program

制御対象

位相補償器

開ループ伝達関数

閉ループ伝達関数

制御前後のステップ応答


%系の個有値の計算

%program

M=0.628;
eta=4.01;
G=3.18;
m=4.735*0.001;
J=1.155*0.0001;
l=0.2465;
xi=1.003*0.0001;
g=9.8;
delta=J*(M+m)+M*m*l*l;
x1=-m*m*g*l*l/delta;
x2=-eta*(J+m*l*l)/delta;
x3=m*l*xi/delta;
x4=m*g*l*(M+m)/delta;
x5=m*l*eta/delta;
x6=-xi*(M+m)/delta;
A=[0,0,1,0;0,0,0,1;0,x1,x2,x3;0,x4,x5,x6]
[v,D]=eig(A)
B=[0;0;(J+m*l*l)*G/delta;-m*l*G/delta]

%Results


A =



         0         0    1.0000         0

         0         0         0    1.0000

         0   -0.0526   -6.3716    0.0005

         0   28.5205   18.4440   -0.2501





v =



    1.0000   -0.0002   -0.0019   -0.0148

         0    0.1885   -0.1830    0.1523

         0   -0.0008    0.0101    0.0957

         0    0.9821    0.9831   -0.9836





D =



         0         0         0         0

         0    5.2094         0         0

         0         0   -5.3711         0

         0         0         0   -6.4600



B =

         0

         0

    5.0528

  -14.6264

極配置法によるfeedbackgainの決定

%Program

s1=-14;
s2=-18;
s3=-12;
s4=-20;
g1=3;
g2=5;
g3=7;
g4=9;
I=[1 0 0 0;0 1 0 0;0 0 1 0; 0 0 0 1];
w1=inv(s1*I-A)*B*g1
w2=inv(s2*I-A)*B*g2
w3=inv(s3*I-A)*B*g3
w4=inv(s4*I-A)*B*g4
K=-[g1 g2 g3 g4]*inv([w1 w2 w3 w4])

%Results


w1 =

    0.1422

   -0.4915

   -1.9909

    6.8813



w2 =

    0.1208

   -0.3892

   -2.1746

    7.0051



w3 =

    0.5253

   -1.9440

   -6.3042

   23.3281



w4 =

    0.1670

   -0.5272

   -3.3392

   10.5450



K =

 -421.9381 -261.1163 -114.7987  -43.5809

%Program:

clf

A1=A-B*K
B1=B;
C1=[0 1 0 0];
D1=0;
impulse(A1,B1,C1,D1)

%Results:


A-BK =

  1.0e+003 *

         0         0    0.0010         0

         0         0         0    0.0010

    2.1320    1.3193    0.5737    0.2202

   -6.1714   -3.7907   -1.6606   -0.6377



%Program

clf
A1=A-B*K
B1=B;
C1=[1 0 0 0];
D1=0;
impulse(A1,B1,C1,D1)

%Results





最適regulatorの構成

Case 1 : diag(1,1,1,1)
Program

A=[0 0 1 0;0 0 0 1;0 -0.0526 -6.3716 0.0005;0 28.5205 18.4440 -0.2501];
B=[0;0;5.0258;-14.6264];
C=[1 0 0 0;0 1 0 0];
D=[0;0];
r=1
x0=[0;0.1;0;0];
Q=[1 0 0 0;0 1 0 0;0 0 1 0;0 0 0 1]

Kp1=lqr(A,B,r,Q)

t=0.0:0.005:5.0;

y=initial(A-B*Kp1,B,[C;-Kp1],[D;0],x0,t);

subplot(3,1,1), plot(t,y(:,1)),grid
xlabel('Time[sec]'), ylabel('z:Amplitude[m]')

subplot(3,1,2), plot(t,y(:,2)),grid
xlabel('Time[sec]'), ylabel(' :Amplitude[rad]')

subplot(3,1,3), plot(t,y(:,3)),grid
xlabel('Time[sec]'), ylabel('u:Input[v]')

%Results:


r =



     1





Q =



     1     0     0     0

     0     1     0     0

     0     0     1     0

     0     0     0     1





Kp1 =



   -1.0000  -12.0136   -3.3915   -2.4127




Case 2 : diag(100,1,1,1)
%Program

A=[0 0 1 0;0 0 0 1;0 -0.0526 -6.3716 0.0005;0 28.5205 18.4440 -0.2501];
B=[0;0;5.0258;-14.6264];
C=[1 0 0 0;0 1 0 0];
D=[0;0];
r=1
x0=[0;0.1;0;0];
Q=[100 0 0 0;0 1 0 0;0 0 1 0;0 0 0 1]

Kp1=lqr(A,B,r,Q)

t=0.0:0.005:5.0;

y=initial(A-B*Kp1,B,[C;-Kp1],[D;0],x0,t);

clf

subplot(3,1,1), plot(t,y(:,1)),grid
xlabel('Time[sec]'), ylabel('z:Amplitude[m]')

subplot(3,1,2), plot(t,y(:,2)),grid
xlabel('Time[sec]'), ylabel(' :Amplitude[rad]')

subplot(3,1,3), plot(t,y(:,3)),grid
xlabel('Time[sec]'), ylabel('u:Input[v]')

%Results


r =



     1





Q =



   100     0     0     0

     0     1     0     0

     0     0     1     0

     0     0     0     1





Kp1 =



  -10.0000  -22.3180   -7.8762   -4.2113




Case 3 : diag(1,100,1,1)
Program

A=[0 0 1 0;0 0 0 1;0 -0.0526 -6.3716 0.0005;0 28.5205 18.4440 -0.2501];
B=[0;0;5.0258;-14.6264];
C=[1 0 0 0;0 1 0 0];
D=[0;0];
r=1
x0=[0;0.1;0;0];
Q=[1 0 0 0;0 100 0 0;0 0 1 0;0 0 0 1]

Kp1=lqr(A,B,r,Q)

t=0.0:0.005:5.0;

y=initial(A-B*Kp1,B,[C;-Kp1],[D;0],x0,t);

subplot(3,1,1), plot(t,y(:,1)),grid
xlabel('Time[sec]'), ylabel('z:Amplitude[m]')

subplot(3,1,2), plot(t,y(:,2)),grid
xlabel('Time[sec]'), ylabel(' :Amplitude[rad]')

subplot(3,1,3), plot(t,y(:,3)),grid
xlabel('Time[sec]'), ylabel('u:Input[v]')

%Results


r =



     1

Q =

     1     0     0     0

     0   100     0     0

     0     0     1     0

     0     0     0     1



Kp1 =

   -1.0000  -16.7944   -3.5912   -2.6635




Case 4 : diag(1,1,100,1)
Program

A=[0 0 1 0;0 0 0 1;0 -0.0526 -6.3716 0.0005;0 28.5205 18.4440 -0.2501];
B=[0;0;5.0258;-14.6264];
C=[1 0 0 0;0 1 0 0];
D=[0;0];
r=1
x0=[0;0.1;0;0];
Q=[1 0 0 0;0 1 0 0;0 0 100 0;0 0 0 1]

Kp1=lqr(A,B,r,Q)

t=0.0:0.005:5.0;

y=initial(A-B*Kp1,B,[C;-Kp1],[D;0],x0,t);

subplot(3,1,1), plot(t,y(:,1)),grid
xlabel('Time[sec]'), ylabel('z:Amplitude[m]')

subplot(3,1,2), plot(t,y(:,2)),grid
xlabel('Time[sec]'), ylabel(' :Amplitude[rad]')

subplot(3,1,3), plot(t,y(:,3)),grid
xlabel('Time[sec]'), ylabel('u:Input[v]')

%Results


r =



     1





Q =



     1     0     0     0

     0     1     0     0

     0     0   100     0

     0     0     0     1





Kp1 =



   -1.0000  -42.5051  -11.6069   -7.8573




Case 5 : diag(1,1,1,100)
Program

A=[0 0 1 0;0 0 0 1;0 -0.0526 -6.3716 0.0005;0 28.5205 18.4440 -0.2501];
B=[0;0;5.0258;-14.6264];
C=[1 0 0 0;0 1 0 0];
D=[0;0];
r=1
x0=[0;0.1;0;0];
Q=[1 0 0 0;0 1 0 0;0 0 1 0;0 0 0 100]

Kp1=lqr(A,B,r,Q)

t=0.0:0.005:5.0;

y=initial(A-B*Kp1,B,[C;-Kp1],[D;0],x0,t);

subplot(3,1,1), plot(t,y(:,1)),grid
xlabel('Time[sec]'), ylabel('z:Amplitude[m]')

subplot(3,1,2), plot(t,y(:,2)),grid
xlabel('Time[sec]'), ylabel(' :Amplitude[rad]')

subplot(3,1,3), plot(t,y(:,3)),grid
xlabel('Time[sec]'), ylabel('u:Input[v]')

%Results


r =



     1





Q =



     1     0     0     0

     0     1     0     0

     0     0     1     0

     0     0     0   100





Kp1 =



   -1.0000  -25.9687   -3.9366  -11.0898






〜前期中間試験〜

問3 program

根軌跡

問4 program

プラント

コントローラ

ボード線図