MAE 280 B
Linear Control Design

LMI LQR Control

Mauricio de Oliveira & Chris Nelson

<<NC` 
<<NCAlgebra`

<<NCSDP`
<<SDPSylvester`

LMI LQR Control

Satellite problem data

m = 100 (* 100 kg *);
r = 300*10^3 (* 300 km *);
R0 = 6.37*10^6 (* 6.37 10^3 km *);
G = 6.673*10^(-11) (* 6.673 N m^2/kg^2 *);
M = 5.98*10^(24) (* 5.98 10^24 kg *);
k = G*M; (* gravitational force constant *)
w = Sqrt[k/((R0+r)^3)];(* angular velocity (rad/s) *)
v = w*(R0+r);(* "ground" velocity (m/s) *)

Linearlized system matrices

A = {{0,0,1,0},{0,0,0,1},{3*w^2,0,0,2*(r+R0)*w},{0,0,-2*w/(r+R0),0}};
Bu = {{0,0},{0,0},{1/m,0},{0,1/(m*r)}};
Bw = {{0,0},{0,0},{1/m,0},{0,1/(m*r)}};

Noise variances

W = 0.1*IdentityMatrix[2];

Scale

T = DiagonalMatrix[{1,r,1,r}]

Similarity transformation

At = T.A.Inverse[T]
But = T.Bu
Bwt = T.Bw

n = Length[At];

Optimal LQR state feedback control (using u2 only)

LQR Solution

m = 1;

Cz = Join[IdentityMatrix[n],ConstantArray[0,{m,n}]];
Dzu = Join[ConstantArray[0,{n,m}],IdentityMatrix[m]];

K = LQRegulatorGains[StateSpaceModel[{At,But[[All,{2}]]}],{Transpose[Cz].Cz,Transpose[Dzu].Dzu}];
Klqr = -K

LMI Solution

SNC[a,b,c,d,bw,q,q1,x,l,z];
F={
{{a**x+b**l+x**tp[a]+tp[b**l],tp[c**x+d**l]},
{c**x+d**l,-1}},
-{{z,tp[bw]},
{bw,x}}
};
vars={z,x,l};
obj={-1,0,0};

NCSDPForm[F,vars,obj]

data={a->At,b->But[[All,{2}]],c->Cz,d->Dzu,bw->Bwt[[All,{2}]]};

{abc,rules}=NCSDP[F,vars,obj,data];

{Y,X,S,iters}=SDPSolve[abc,rules,GapTol->10.^(-18)];

Check the inequalities

{Z,X,L}=Y;

LMI1=ArrayFlatten[{{At.X+But[[All,{2}]].L,ConstantArray[0,{n,m+n}]},{Cz.X+Dzu.L,-IdentityMatrix[m+n]}}];
LMI1=.5*(LMI1+Transpose[LMI1]);
Eigenvalues[LMI1]

LMI2=ArrayFlatten[{{Z,Transpose[Bwt[[All,{2}]]]},{Bwt[[All,{2}]],X}}];
Eigenvalues[LMI2]

Construct K and check closed loop stability

K = L.Inverse[X]

Acl=At+But[[;;,2;;2]].K;
Eigenvalues[Acl]

Klqr

Optimal LQR state feedback control (using u1 and u2)

LQR Solution

m = 2;

Cz = Join[IdentityMatrix[n],ConstantArray[0,{m,n}]];
Dzu = Join[ConstantArray[0,{n,m}],IdentityMatrix[m]];

K = LQRegulatorGains[StateSpaceModel[{At,But}],{Transpose[Cz].Cz,Transpose[Dzu].Dzu}];
Klqr = -K

LMI Solution

SNC[a,b,c,d,bw,x,l,z];
F={
{{a**x+b**l+x**tp[a]+tp[b**l],tp[c**x+d**l]},
{c**x+d**l,-1}},
-{{z,tp[bw]},
{bw,x}}
};
vars={z,x,l};
obj={-1,0,0};

NCSDPForm[F,vars,obj]

data={a->At,b->But,c->Cz,d->Dzu,bw->Bwt};

{abc,rules}=NCSDP[F,vars,obj,data,UserRules->{tp[q]->q}];

{Y,X,S,iters}=SDPSolve[abc,rules];

Check the inequalities

{Z,X,L}=Y;

LMI1=At.X+But.L;
LMI1=.5*(LMI1+Transpose[LMI1]);
Eigenvalues[LMI1]

LMI2=ArrayFlatten[{{Z,Transpose[Bwt]},{Bwt,X}}];
Eigenvalues[LMI2]

Construct K and check closed loop stability

K = L.Inverse[X]

Klqr

Acl=At+But.K;
Eigenvalues[Acl]