ŘMS 1: samostatná práce


Simulační model

Jedná se o matematický model inverzního kyvadla, odvozený metodou Lagrangeových rovnic II. druhu.

Funkce f(x,u) pravých stran nelineárního systému x’=f(x,u), y=g(x,u) (v samostatném souboru):

function dx = oderightfun(t,x)
%RMS I 2007

dx=zeros(4,1);  %Výstupem je sloupcovy vektor
u=0;   %Vstupní síla
m=1;  %kg
M=1;  %kg
l=0.1; %m
g=9.81; %m/s/s
b=0;

x1=x(1);
x2=x(2);
x3=x(3);
x4=x(4);

dx(1)=x2;
dx(2)=(m*l*sin(x3)*x4^2+u-b*x2-m*g*sin(x3)*cos(x3))/(m*sin(x3)^2+M);
dx(3)=x4;
dx(4)=(g*sin(x3)*(m+M)-sin(x3)*cos(x3)*x4^2-u*cos(x3)+b*cos(x3)*x2)/(l*(m*sin(x3)^2+M));

Simulační ověření nelineárního modelu:

%Simulace nelinearniho modelu systemu - vyzaduje funkci-soubor oderightfun.m
%Integrace s nulovymi p.p., vykresleni prubehu stavu
ode45('oderightfun',[0 10],[0 0 0 0])
%Integrace s nenulovymi p.p., ulozeni prubehu casu a stavu do vektoru t,x
[t,x]= ode45('oderightfun',[0 10],[0 0 0.001 0]);
%Vykresleni vysledku
plot(t,x)
grid
zoom on
legend('x1','x2','x3','x4')

Linearizovaný model pro rovnovážný stav x=[0 0 0 0]:

%Definice matic linearniho systemu - z linearizace pomoci Taylorova rozvoje:
u=0;   %Vstupní síla
m=1;  %kg
M=1;  %kg
l=0.1; %m
g=9.81; %m/s/s
b=0.01;
u=0;
A= [0 1 0 0; 0 0 -m*g/M 0; 0 0 0 1; 0 0 g*(m+M)/(M*l) 0];
B=[0 1/M 0 -1/(l*M)]';


Analýza linearizovaného modelu

Stabilita lineárního systému

eig(A)
ans =
                  0
                  0
  14.00714103591450
 -14.00714103591450

Póly nejsou v levé polorovině, systém je nestabilní.

Analýza řiditelnosti

rank(ctrb(A,B))
ans =
     4

Hodnost matice řiditelnosti odpovídá řádu systému, systém je řiditelný.

Vhodná volba výstupu pro zajištění pozorovatelnosti:

%Jako vystup volíme uhel natoceni kyvadla
C=[0 0 1 0];
D=[0];
rank(obsv(A,C))
ans =
     2
%      S takovym vystupem system neni pozorovatelný

%Jako vystup volíme miru odpovidajici poloze koncoveho bodu kyvadla
C=[1 0 1 0];
D=[0];
rank(obsv(A,C))
ans =
     4
%System pozorovatelný je.

Určení vlastních frekvencí a poměrného tlumení systému:

>> damp(A)                                            
 Eigenvalue      Damping     Freq. (rad/s)                                             
  1.40e+001    -1.00e+000      1.40e+001    
  0.00e+000    -1.00e+000      0.00e+000    
  0.00e+000    -1.00e+000      0.00e+000    
 -1.40e+001     1.00e+000      1.40e+001

Zjištění odezvy na jednotkový skok a Diracův impuls:

figure; step(A,B,C,D);
figure; impulse(A,B,C,D);

Systém je nestabilní, proto graf odezvy na Diracův impuls i jednotkový skok diverguje.

Nyquistův diagram

figure; nyquist(A,B,C,D);

Bode diagram

figure; bode(A,B,C,D);


Návrh řízení a simulační ověření


P-regulátor a Geometrické místo kořenů

Zapojíme zápornou zpětnou vazbu (viz help rlocus):

                      +-----+
          ---->O----->| SYS |----+---->
              -|      +-----+    |
               |                 |
               |       +---+     |
               +-------| K |<----+
                       +---+

Vyšetříme polohu kořenů systému při změně K=[0:inf]:

rlocus(A,B,C,D);

Ze získaného grafu vyplývá (póly se nepodaří přemístit do levé poloroviny), že vyšetřovaný systém není stabilizovatelný P regulátorem.


Stavová zpětná vazba navržená metodou umisťování pólů

Volíme žádanou polohu pólů

P=[-1 -1 -1 -1];

Vypočteme příslušné zesílení stavové zpětné vazby

K=acker(A,B,P);

Kontrolujeme stabilitu řízeného lineárního systému

eig(A-B*K);
ans =
 -1.00054882365801                    
 -0.99999984798141 + 0.00054867159806i
 -0.99999984798141 - 0.00054867159806i
 -0.99945148037917  

Linearizovaný systém je navrženou stavovou zpětnou vazbou stabilizován, póly jsou umístěny s malou tolerancí podle požadavku.

Navržené řízení kontrolujeme pomocí přechodových charakteristik odezvy na jednotkový skok a Diracův impuls pro nový vstup w (u=K*x+w) a frekvenční charakteristiky (Nyquistův diagram):

figure; step(A-B*K,B,C,D);
figure; impulse(A-B*K,B,C,D);
figure; nyquist(A-B*K,B,C,D);

Navržený regulátor zapojíme do obvodu s modelem původního nelineárního systému a kontrolujeme simulací v časové oblasti:

Vytvoříme funkci pravých stran se stavovou zpětnou vazbou

function dx = oderightfun_ctrb(t,x)
%RMS I 2007

dx=zeros(4,1);  %Sloupcovy vektor

%Stavová zpětná vazba 
K=[-0.01019367991845  -0.04077471967380 -20.22101936799184  -0.40407747196738];
u=-K*x;   %Vstupní síla

m=1;  %kg
M=1;  %kg
l=0.1; %m
g=9.81; %m/s/s
b=0;

x1=x(1);
x2=x(2);
x3=x(3);
x4=x(4);

dx(1)=x2;
dx(2)=(m*l*sin(x3)*x4^2+u-b*x2-m*g*sin(x3)*cos(x3))/(m*sin(x3)^2+M);
dx(3)=x4;
dx(4)=(g*sin(x3)*(m+M)-sin(x3)*cos(x3)*x4^2-u*cos(x3)+b*cos(x3)*x2)/(l*(m*sin(x3)^2+M));

Simulační ověření řízeného nelineárního modelu:

ode45('oderightfun_ctrb',[0 10],[0.01 0 0 0]);


Stavová zpětná vazba navržená LQR (pro nekonečný čas ustálení)

Volíme pozitivně definitní váhové matice Q,R kriteria J = Integral {x’Qx + u’Ru} dt, např.

Q=diag([1e3,1e3,1e3,1e3])
R=1

Vypočteme zesílení stavové zpětné vazby:

K=lqr(A,B,Q,R)
K =
  1.0e+002 *
  -0.31622776601686  -0.46827275425043  -2.04633899588985  -0.37090972240201

a dále postupujeme jako v případě návrhu pomocí Ackermanovy formule.


Stavový pozorovatel

Navrhujeme stavového pozorovatele, který na základě měřených výstupů řízeného systému y a počátečních podmínek x0 rekonstruuje odhad stavového vektoru xst.

%xst' = (A-L*C)xst + B*u + L*y

Koeficienty zesílení zpětné vazby pozorovatele L vypočteme (pro požadované umístění pólů -1*[11 12 13 14])

 L = acker(A',C',-1*[11 12 13 14])'

Připravíme si objekt našeho stavového popisu ve formátu Control Toolboxu

sys_pp = ss(A,B,C,D)

Příkazem reg vytvoříme objekt výstupního regulátoru, který obsahuje stavový pozorovatel a stavový regulátor (K*xst), vypočtený v předchozích krocích.

 sys_reg = reg(sys_pp,K,L) 

Tento výstupní regulátor zapojíme do KLADNÉ zpětné vazby s původním systémem:

 sys_controlled = feedback(sys_pp,sys_reg,+1) 

Celý model má nyní čtyři vlastní stavy a čtyři stavy pozorovatele. Výpočet odezvy pro soustavu s nenulovými počátečními podmínkami [0 0 0.1 0] vlastního modelu a [0 0 0 0] pozorovatele:

 
[Y,T,X]=initial(sys_controlled, [0 0 0.1 0 0 0 0 0]);
plot(T, X(:,3), T, X(:,7))

Stránky odboru mechaniky a mechatroniky FS ČVUT v Praze