第3章基于LMI的输入受限控制算法设计

3.1控制输入受限下的LMI控制算法设计
在实际控制系统的设计中，通常面临控制输入受限的问题。采用LMI方法，可实现控制输入受限下的控制算法设计［1~4］。
3.1.1系统描述
考虑状态方程
x·=Ax+Bu(3.1)
其中x=［x1x2］T，u为控制输入。
控制器设计为
u=Kx(3.2)
其中，K=［k1k2］。
控制目标为通过设计LMI求解K，实现x→0，|u|≤umax。
3.1.2控制器的设计与分析
设计Lyapunov函数如下： 

V=xTPx

其中，P>0，P=PT。通过P的设计可有效地调节x的收敛效果，并有利于LMI的求解。
则
V·=x·TPx+xTPx·
=（Ax+Bu）TPx+xTP（Ax+Bu）
=（Ax+BKx）TPx+xTP（Ax+BKx）
=xT（A+BK）TPx+xTP（A+BK）x
=xTQT1x+xTQ1x=xTQx
其中，Q1=P（A+BK），Q=QT1+Q1。
为了实现指数收敛，即V·≤－αV，取
αV+V·=αxTPx+xTQx=xT（αP+Q）x
取αP+Q<0，α>0，则αV+V·≤0，即V·≤－αV，由V·≤－αV，可得解为
V（t）≤V（0）exp（－αt）≤V（0）
如果t→∞, 则V（t）→0，从而x→0且指数收敛。
3.1.3LMI的设计与求解
由于V（0）=xT0Px0，如果存在正定对称阵P和ω->0，使得xT0Px0≤ω-成立，则可保证V（0）≤ω-，从而保证V（t）≤ω-。
取KTK≤ω-－1u2maxP，由u=Kx可得
u2=（Kx）TKx=xTKTKx≤xTω-－1u2maxPx=ω-－1u2maxV≤u2max
则
|u|≤umax
通过上述分析，构造两个LMI为
αP+Q<0(3.3)

KTK－ω-－1u2maxP≤0(3.4)
不等式(3.3)中，Q中含有P和K，式(3.4)中也含有P和K，故不能独立Q存在，将式(3.3)中的Q展开如下： 
αP+PA+PBK+ATP+KTBTP<0
左右同乘以P－1，可得
αP－1+AP－1+BKP－1+P－1AT+P－1KTBT<0(3.5)
不等式(3.4)中含有非线性项，必须转化为线性矩阵不等式才能求解。取k0=ω-－1u2max，则不等式(3.4)变为KTK≤k0P。根据Schur补定理，式(3.4)变换为
k0PKT
K1≥0
左右同乘以P－10
01，可得
k0P－1P－1KT
KP－11≥0(3.6)
考虑式(3.5)和式(3.6),令F=KP－1和 N=P－1,则P－1KT=FT，由式(3.5)和式(3.6)可得第1个和第2个LMI为
αN+AN+BF+NAT+FTBT<0(3.7)

k0NFT
F1≥0(3.8)
根据P的定义，可设计第3个LMI为
P>0，P=PT(3.9)
要满足xT0Px0≤ω-，根据Schur补定理，可设计第4个LMI为
ω-xT0
x0N≥0(3.10)
通过上面4个LMI，即式(3.7)至式(3.10)，通过设计合适的umax和α值，可求得有效的K。
3.1.4仿真实例
实际模型为
θ¨=－25θ·+133u（t）
考虑模型式(3.1)，取x1=θ，x2=θ·，则对应于式x·=Ax+Bu，有A=01
0－bJ，B=0
1J。取J=1133，b=25133，初始状态值为x（0）=［10］。
取ω-=0.10，α=10，umax=1.0，采用LMI程序chap3_1LMI.m，求解LMI式(3.7)至式(3.10)，MATLAB运行后显示有可行解，解为K=［－0.9528－0.0381］。控制律采用式(3.2)，将求得的K代入控制器程序chap3_1ctrl.m，仿真结果如图3.1和图3.2所示。为了保证有可行解，可取较大的umax值，并取α为较小的值。



图3.1状态响应





图3.2控制输入信号


仿真程序： 
(1) LMI不等式求K程序： chap3_1LMI.m

clear all;

close all;



J=1/133;b=25/133;

A=［0 1;0 -b/J］;

B=［0 1/J］';



F=sdpvar(1,2);

P=sdpvar(2,2,'symmetric');

N=sdpvar(2,2,'symmetric');



umax=1.0;

alfa=10;w_bar=0.10;

x0=［1 0］';



%First LMI

L1=set((alfa*N+A*N+B*F+N*A'+F'*B')<0);



%Second LMI

k0=umax^2/w_bar;

M1=［k0*N F';F 1］; 

L2=set(M1>0);



%Third LMI

L3=set(N>0);



%Fourth LMI

M2=［w_bar x0';x0 N］;

L4=set(M2>0);





L=L1+L2+L3+L4;

solvesdp(L);



F=double(F);

N=double(N);



P=inv(N)

K=F*P

(2)  Simulink主程序： chap3_1sim.mdl





(3) 被控对象S函数： chap3_1plant.m

function ［sys,x0,str,ts］ = spacemodel(t,x,u,flag)

switch flag,

case 0,

［sys,x0,str,ts］=mdlInitializeSizes;

case 1,

sys=mdlDerivatives(t,x,u);

case 3,

sys=mdlOutputs(t,x,u);

case {2,4,9}

sys=［］;

otherwise

error(［'Unhandled flag = ',num2str(flag)］);

end

function ［sys,x0,str,ts］=mdlInitializeSizes

sizes = simsizes;

sizes.NumContStates= 2;

sizes.NumDiscStates  = 0;

sizes.NumOutputs     = 2;

sizes.NumInputs      = 1;

sizes.DirFeedthrough = 0;

sizes.NumSampleTimes = 0;

sys = simsizes(sizes);

x0 =［1 0］;

str = ［］;

ts = ［］;

function sys=mdlDerivatives(t,x,u)

A=［0 1;0 -25］;

B=［0 133］';



ut=u(1);



dx=A*x+B*ut;



sys(1)=dx(1);

sys(2)=dx(2);

function sys=mdlOutputs(t,x,u)

sys(1)=x(1);

sys(2)=x(2);

(4) 控制器S函数： chap3_1ctrl.m

function ［sys,x0,str,ts］ = spacemodel(t,x,u,flag)

switch flag,

case 0,

［sys,x0,str,ts］=mdlInitializeSizes;

case 3,

sys=mdlOutputs(t,x,u);

case {2,4,9}

sys=［］;

otherwise

error(［'Unhandled flag = ',num2str(flag)］);

end

function ［sys,x0,str,ts］=mdlInitializeSizes

sizes = simsizes;

sizes.NumContStates= 0;

sizes.NumDiscStates  = 0;

sizes.NumOutputs     = 1;

sizes.NumInputs      = 2;

sizes.DirFeedthrough = 1;

sizes.NumSampleTimes = 1;

sys = simsizes(sizes);

x0  = ［］;

str = ［］;

ts  = ［0 0］;

function sys=mdlOutputs(t,x,u)

x1=u(1);

x2=u(2);

X=［x1 x2］';



K=［ -0.9528   -0.0381］;



ut=K*X;

sys(1)=ut;

(5) 作图程序： chap3_1plot.m

close all;



figure(1);

subplot(211);

plot(t,x(:,1),'r','linewidth',2);

xlabel('time(s)');ylabel('x1 response');

subplot(212);

plot(t,x(:,2),'b','linewidth',2);

xlabel('time(s)');ylabel('x2 response');



figure(2);

plot(t,ut(:,1),'r','linewidth',2);

xlabel('time(s)');ylabel('ut');

3.2控制输入受限下位置跟踪LMI控制算法
3.2.1系统描述

考虑电机负载模型如下： 
Jθ¨=－bθ·+u（t）(3.11)
其中，θ为角度，J为转动惯量，b为黏性系数，u为控制输入，|u（t）|≤umax。
式(3.11)可写为
θ¨=－bJθ·+1Ju（t）
取角度指令为θd，则角度跟踪误差为x1=θ－θd，角速度跟踪误差为x2=θ·－θ·d，则控制目标在控制输入受限条件下，实现角度和角速度跟踪，即t→∞时，x1→0，x2→0，|u（t）|≤umax。
由于
x·2=－bJθ·+1Ju－θ¨d=－bJ（x2+θ·d）+1Ju－θ¨d
=－bJx2+1Ju－θ¨d－bJθ·d
=－bJx2+1J（u－Jθ¨d－bθ·d）
取τ=u－Jθ¨d－bθ·d，即
u=τ+Jθ¨d+bθ·d(3.12)
则|u（t）|≤umax转化为|τ+Jθ¨d+bθ·d|≤umax，即－umax≤τ+Jθ¨d+bθ·d≤umax，从而
τ≤umax－Jθ¨d－bθ·d≤umax+J|θ¨d|max+b|θ·d|max
即
τmax=umax+J|θ¨d|max+b|θ·d|max
从而由umax可得到τmax。
由x·2=－bJx2+1Jτ，可得
x·1=x2
x·2=－bJx2+1Jτ
则误差状态方程为
x·=Ax+Bτ(3.13)
其中，x=［x1x2］T，A=00
0－bJ，B=0
1J。
控制器设计为
τ=Kx(3.14)
其中，K=［k1k2］。
控制目标转化为通过设计LMI，实现x→0，|τ|≤τmax。
针对模型式(3.13)进行控制器的设计、收敛性分析及LMI的设计，与3.1节相同。
3.2.2仿真实例
实际模型为
θ¨=－25θ·+133u（t）
被控对象角度和角速度的初始值为［1.00］T，取角度指令为θd=sint，则θ·d=cost，角度跟踪误差为x1（0）=θ（0）－θd（0）=1.0，角速度跟踪误差为x2（0）=θ·（0）－θ·d（0）=－1.0，x（0）=［1－1］。
取ω-=0.10，α=10，umax=1.0，采用LMI程序chap3_2LMI.m，求解LMI式(3.7)至式(3.10)，MATLAB运行后显示有可行解，解为K=［－0.9870－0.0293］。控制律采用式(3.12)，将求得的K代入控制器程序chap3_2ctrl.m，仿真结果如图3.3和图3.4所示。为了保证有可行解，可取umax为较大的值，并取α为较小的值。



图3.3角度和角速度跟踪





图3.4控制输入信号




仿真程序： 
(1) LMI不等式求K程序： chap3_2LMI.m

clear all;

close all;

 

J=1/133;

b=25/133;

 

A=［0 0;0 -b/J］;

B=［0 1/J］';

 

K=sdpvar(1,2);

 

M=sdpvar(3,3);

F=sdpvar(1,2);

 

P=sdpvar(2,2,'symmetric');

N=sdpvar(2,2,'symmetric');

 

umax=1.0;

tol_max=umax+J+b;

alfa=10;w_bar=1.0;

 

%First LMI

x0=［1 -1］';

L1=set((alfa*N+A*N+B*F+N*A'+F'*B')<0);

 

%Second LMI

k0=tol_max^2/w_bar;

M=［k0*N F';F 1］; 

L2=set(M>0);

 

%Third LMI

L3=set(N>0);

 

%Fourth LMI

M1=［w_bar x0';x0 N］;

L4=set(M1>0);

 

 

L=L1+L2+L3+L4;

solvesdp(L);

 

F=double(F);

N=double(N);

 

P=inv(N)

K=F*P

(2)  Simulink主程序： chap3_2sim.mdl





(3) 被控对象S函数： chap3_2plant.m

function ［sys,x0,str,ts］ = spacemodel(t,x,u,flag)

switch flag,

case 0,

［sys,x0,str,ts］=mdlInitializeSizes;

case 1,

sys=mdlDerivatives(t,x,u);

case 3,

sys=mdlOutputs(t,x,u);

case {2,4,9}

sys=［］;

otherwise

error(［'Unhandled flag = ',num2str(flag)］);

end

function ［sys,x0,str,ts］=mdlInitializeSizes

sizes = simsizes;

sizes.NumContStates= 2;

sizes.NumDiscStates  = 0;

sizes.NumOutputs     = 2;

sizes.NumInputs      = 1;

sizes.DirFeedthrough = 0;

sizes.NumSampleTimes = 0;

sys = simsizes(sizes);

x0 =［1 0］;

str = ［］;

ts = ［］;

function sys=mdlDerivatives(t,x,u)

A=［0 1;0 -25］;

B=［0 133］';



ut=u(1);

 

dx=A*x+B*ut;

 

sys(1)=dx(1);

sys(2)=dx(2);

function sys=mdlOutputs(t,x,u)

sys(1)=x(1);

sys(2)=x(2);

(4) 控制器S函数： chap3_2ctrl.m

function ［sys,x0,str,ts］ = spacemodel(t,x,u,flag)

switch flag,

case 0,

［sys,x0,str,ts］=mdlInitializeSizes;

case 3,

sys=mdlOutputs(t,x,u);

case {2,4,9}

sys=［］;

otherwise

error(［'Unhandled flag = ',num2str(flag)］);

end

function ［sys,x0,str,ts］=mdlInitializeSizes

sizes = simsizes;

sizes.NumContStates= 0;

sizes.NumDiscStates  = 0;

sizes.NumOutputs     = 1;

sizes.NumInputs      = 2;

sizes.DirFeedthrough = 1;

sizes.NumSampleTimes = 1;

sys = simsizes(sizes);

x0  = ［］;

str = ［］;

ts  = ［0 0］;

function sys=mdlOutputs(t,x,u)

x1=u(1);

x2=u(2);

 

xd=sin(t);

dxd=cos(t);

ddxd=-sin(t);

e=x1-xd;

de=x2-dxd;

 

K=［-0.9870   -0.0293］;

 

tol=K*［e de］';

 

ut=tol+1/133*ddxd+25/133*dxd;

 

sys(1)=ut;

(5) 作图程序： chap3_2plot.m

close all;

 

figure(1);

subplot(211);

plot(t,sin(t),'r',t,x(:,1),'b','linewidth',2);

xlabel('time(s)');ylabel('Angle tracking');

subplot(212);

plot(t,cos(t),'r',t,x(:,2),'b','linewidth',2);

xlabel('time(s)');ylabel('Angle speed tracking');

 

figure(2);

plot(t,ut(:,1),'r','linewidth',2);

xlabel('time(s)');ylabel('ut');

3.3控制输入受限下带扰动控制系统LMI控制算法设计
3.3.1系统描述

考虑状态方程
x·=Ax+B（u+d）(3.15)
其中，x=［x1x2］T，u为控制输入，|u|≤umax，d（t）为加在控制输入上的积分有界扰动。
控制器设计为
u=Kx(3.16)
其中，K=［k1k2］。
控制目标为在满足|u|≤umax条件下，通过设计LMI求解K，实现x→0。
3.3.2基于H∞指标控制器的设计与分析
设计Lyapunov函数如下： 
V=xTPx
其中，P>0，P=PT。
通过P的设计，可有效地调节x的收敛效果，并有利于LMI的求解。则
V·=x·TPx+xTPx·
=（Ax+Bu+Bd）TPx+xTP（Ax+Bu+Bd）
=（Ax+BKx+Bd）TPx+xTP（Ax+BKx+Bd）
=xT（A+BK）TPx+xTP（A+BK）x+（Bd）TPx+xTP（Bd）
=xTQT1x+xTQ1x+（BTPx+xTPB）d
其中，Q1=P（A+BK），Q=QT1+Q1。
令η=［xTd］T，则
η=x1
x2
d=x
d，ηT=［xTd］=［x1x2d］
从而
V·=xTQx+ηT0PB
BTP0η=ηTQPB
BTP0η
其中，
xTQx=［xTd］Q0
00x
d=ηTQ0
00η

（BTPx+xTPB）d=［dBTPxTPB］x
d
=［xTd］0PB
BTP0x
d=ηT0PB
BTP0η
输出为Z=Cx，H∞指标取
∫t0ZTZdt<∫t0γ2d2（t）dt+V（0）(3.17)
其中，γ>0，C=10
01。
由于
ZTZ－γ2d2=xTCTCx－γ2d2

ηTCTC0
0－γ2η=
［xTd］CTC0
0－γ2x
d

=［xTCTC－γ2d］x
d=xTCTCx－γ2d2
则
ZTZ－γ2d2=ηTCTC0
0－γ2η
从而
V·+ZTZ－γ2d2=ηTQ+CTCPB
（PB）T－γ2η
取
θ=Q+CTCPB
（PB）T－γ2<0(3.18)
则
V·+ZTZ－γ2d2≤0
对上式积分，可得
V+∫t0ZTZdt≤∫t0γ2d2dt+V（0）
假设d为递减的扰动信号，即积分有界扰动［4］，取
∫∞0d2dt≤γ－2vmax
由于∫t0ZTZdt≥0，则
V（t）≤ω-
其中，vmax+V（0）≤ω-。
由V（t）≤ω-可得
Pmin‖x‖2≤xTPx≤ω-
则收敛结果为
‖x‖2≤vmax+V（0）Pmin
3.3.3闭环系统LMI的设计
由式(3.18)展开，可得
PA+PBK+ATP+KTBTPPB
（PB）T－γ2<0
左右同乘以P－10
01，可得
AP－1+BKP－1+P－1AT+P－1KTBTB
BT－γ2<0(3.19)
令F=KP－1，N=P－1，则P－1KT=FT，由式(3.19)，可得第1个LMI为
AN+BF+NAT+FTBTB
BT－γ2<0(3.20)
根据P的定义，可设计第2个LMI为
N>0，N=NT(3.21)
3.3.4控制输入受限下LMI的设计
由于V（0）=xT0Px0，如果存在正定对称阵P和ω->0，使得xT0Px0≤ω-成立，则可保证V（0）≤ω-，从而保证V（t）≤ω-。
取KTK≤ω-－1u2maxP，由u=Kx可得
u2=（Kx）TKx=xTKTKx≤xTω-－1u2maxPx=ω-－1u2maxV≤u2max
则
|u|≤umax
通过上述分析，构造LMI如下： 
KTK－ω-－1u2maxP≤0(3.22)
不等式(3.22)中含有非线性项，必须转化为线性矩阵不等式才能求解。取k0=ω-－1u2max，则不等式(3.22)变为KTK≤k0P。根据Schur补定理，式(3.22)变换为
k0PKT
K1≥0
左右同乘以P－10
01，可得
k0P－1P－1KT
KP－11≥0(3.23)
令F=KP－1和N=P－1,则P－1KT=FT，由式(3.23)，可得第3个LMI为
k0NFT
F1≥0(3.24)
要满足xT0Px0≤ω-，根据Schur补定理，可将其设计为第4个LMI为
ω-xT0
x0N≥0(3.25)
通过上面4个LMI，设计合适的umax，可求得有效的K。
3.3.5仿真实例
实际模型为
θ¨=－25θ·+133［u（t）+d（t）］
考虑模型式(3.1)，取x1=θ，x2=θ·，则对应于式x·=Ax+B（u+d），有A=00
0－bJ，B=0
1J。取J=1133，b=25133，d=0.1e－5t，初始状态值为x（0）=0.10。
取γ=3.0，ω-=0.10，umax=0.50，采用LMI程序chap3_3LMI.m，求解LMI式(3.20)、式(3.21)、式(3.24)和式(3.25)，MATLAB运行后显示有可行解，解为K=－0.00990.1809。控制律采用式(3.16)，将求得的K代入控制器程序chap3_3ctrl.m，仿真结果如图3.5和图3.6所示。



图3.5状态响应





图3.6控制输入信号


仿真程序： 
(1) LMI不等式求K程序： chap3_3LMI.m

clear all;

close all;

 

J=1/133;b=25/133;

A=［0 1;0 -b/J］;

B=［0 1/J］';

 

F=sdpvar(1,2);

P=sdpvar(2,2,'symmetric');

N=sdpvar(2,2,'symmetric');

 

%First LMI

gama=3.0;

M=［A*N+B*F+N*A'+F'*B' B;B' -gama^2］;

L1=set(M<0);

 

%Second LMI

L2=set(N>0);

 

umax=0.50;

w_bar=0.10;

x0=［1 0］';

%Third LMI

k0=umax^2/w_bar;

M1=［k0*N F';F 1］; 

L3=set(M1>0);

%Fourth LMI

M2=［w_bar x0';x0 N］;

L4=set(M2>0);

 

L=L1+L2+L3+L4;

solvesdp(L);

 

F=double(F);

N=double(N);

 

P=inv(N);

K=F*P

(2) Simulink主程序： chap3_3sim.mdl




(3) 被控对象S函数： chap3_3plant.m

function ［sys,x0,str,ts］ = spacemodel(t,x,u,flag)

switch flag,

case 0,

［sys,x0,str,ts］=mdlInitializeSizes;

case 1,

sys=mdlDerivatives(t,x,u);

case 3,

sys=mdlOutputs(t,x,u);

case {2,4,9}

sys=［］;

otherwise

error(［'Unhandled flag = ',num2str(flag)］);

end

function ［sys,x0,str,ts］=mdlInitializeSizes

sizes = simsizes;

sizes.NumContStates= 2;

sizes.NumDiscStates  = 0;

sizes.NumOutputs     = 2;

sizes.NumInputs      = 1;

sizes.DirFeedthrough = 0;

sizes.NumSampleTimes = 0;

sys = simsizes(sizes);

x0 =［0.1 0］;

str = ［］;

ts = ［］;

function sys=mdlDerivatives(t,x,u)

A=［0 1;0 -25］;

B=［0 133］';



ut=u(1);

dt=0.10*exp(-5*t);



dx=A*x+B*(ut+dt);



sys(1)=dx(1);

sys(2)=dx(2);

function sys=mdlOutputs(t,x,u)

sys(1)=x(1);

sys(2)=x(2);

(4) 控制器S函数： chap3_3ctrl.m

function ［sys,x0,str,ts］ = spacemodel(t,x,u,flag)

switch flag,

case 0,

［sys,x0,str,ts］=mdlInitializeSizes;

case 3,

sys=mdlOutputs(t,x,u);

case {2,4,9}

sys=［］;

otherwise

error(［'Unhandled flag = ',num2str(flag)］);

end

function ［sys,x0,str,ts］=mdlInitializeSizes

sizes = simsizes;

sizes.NumContStates= 0;

sizes.NumDiscStates  = 0;

sizes.NumOutputs     = 1;

sizes.NumInputs      = 2;

sizes.DirFeedthrough = 1;

sizes.NumSampleTimes = 1;

sys = simsizes(sizes);

x0  = ［］;

str = ［］;

ts  = ［0 0］;

function sys=mdlOutputs(t,x,u)

x1=u(1);

x2=u(2);

X=［x1 x2］';



K =［-0.0099    0.1809］;

ut=K*X;

sys(1)=ut;

(5) 作图程序： chap3_3plot.m

close all;



figure(1);

subplot(211);

plot(t,x(:,1),'r','linewidth',2);

xlabel('time(s)');ylabel('x1 response');

subplot(212);

plot(t,x(:,2),'b','linewidth',2);

xlabel('time(s)');ylabel('x2 response');



figure(2);

plot(t,ut(:,1),'r','linewidth',2);

xlabel('time(s)');ylabel('ut');

3.4控制输入受限下带扰动控制系统LMI跟踪控制算法设计
3.4.1系统描述

考虑电机负载模型
Jθ¨=－bθ·+u（t）+d（t）(3.26)
其中，θ为角度，J为转动惯量，b为黏性系数，u为控制输入，d（t）为加在控制输入上的积分有界扰动。
式(3.26)可写为
θ¨=－bJθ·+1J［u（t）+d（t）］
取角度指令为θd，则角度跟踪误差为x1=θ－θd，角速度跟踪误差为x2=θ·－θ·d，则控制目标为角度和角速度的跟踪，即t→∞时，x1→0，x2→0。
由于
x·2=－bJθ·+1J（u+d）－θ¨d=－bJ（x2+θ·d）+1J（u+d）－θ¨d
=－bJx2+1J（u+d）－θ¨d－bJθ·d
=－bJx2+1J（u+d－Jθ¨d－bθ·d）
取τ=u－Jθ¨d－bθ·d，即u=τ+Jθ¨d+bθ·d。由x·2=－bJx2+1J（τ+d），可得
x·1=x2

x·2=－bJx2+1J（τ+d）
则误差状态方程为
x·=Ax+B（τ+d）(3.27)
其中，x=［x1x2］T，A=01
0－bJ，B=0
1J。
控制器设计为
τ=Kx(3.28)
其中，K=［k1k2］。
控制目标转化为通过设计LMI求K，实现t→∞时，x→0。
针对模型式(3.27)进行控制器的设计、收敛性分析及LMI的设计，与3.3节“带扰动的控制系统LMI控制算法设计”相同。
3.4.2仿真实例
实际模型为
θ¨=－25θ·+133［u（t）+d（t）］
取x1=θ，x2=θ·，则可得式x·=Ax+B（τ+d），取d=e－5t，初始状态值为［1.00］T。取角度指令为θd=sint，则θ·d=cost，角度跟踪误差为x1（0）=θ（0）－θd（0）=1.0，角速度跟踪误差为x2（0）=θ·（0）－θ·d（0）=－1.0，x（0）=［1－1］。
取γ=3.0，采用LMI程序chap3_4LMI.m，求解LMI式，求解LMI式(3.20)、式(3.21)、式(3.24)和式(3.25)，MATLAB运行后显示有可行解，解为K=［－0.00990.1809］。控制律采用式(3.28)和u=τ+Jθ¨d+bθ·d，将求得的K代入控制器程序chap3_4ctrl.m中，仿真结果如图3.7和图3.8所示。



图3.7角度及角速度跟踪




图3.8控制输入信号


仿真程序： 
(1) LMI不等式求K程序： chap3_4LMI.m

clear all;

close all;

 

J=1/133;b=25/133;

A=［0 1;0 -b/J］;

B=［0 1/J］';

 

F=sdpvar(1,2);

P=sdpvar(2,2,'symmetric');

N=sdpvar(2,2,'symmetric');

 

%First LMI

gama=3.0;

M=［A*N+B*F+N*A'+F'*B' B;B' -gama^2］;

L1=set(M<0);

 

%Second LMI

L2=set(N>0);

 

umax=0.50;

w_bar=0.10;

x0=［1 0］';

%Third LMI

k0=umax^2/w_bar;

M1=［k0*N F';F 1］; 

L3=set(M1>0);

%Fourth LMI

M2=［w_bar x0';x0 N］;

L4=set(M2>0);

 

L=L1+L2+L3+L4;

solvesdp(L);

 

F=double(F);

N=double(N);

 

P=inv(N);

K=F*P

(2)  Simulink主程序： chap3_4sim.mdl




(3) 被控对象S函数： chap3_4plant.m

function ［sys,x0,str,ts］ = spacemodel(t,x,u,flag)

switch flag,

case 0,

［sys,x0,str,ts］=mdlInitializeSizes;

case 1,

sys=mdlDerivatives(t,x,u);

case 3,

sys=mdlOutputs(t,x,u);

case {2,4,9}

sys=［］;

otherwise

error(［'Unhandled flag = ',num2str(flag)］);

end

function ［sys,x0,str,ts］=mdlInitializeSizes

sizes = simsizes;

sizes.NumContStates= 2;

sizes.NumDiscStates  = 0;

sizes.NumOutputs     = 2;

sizes.NumInputs      = 1;

sizes.DirFeedthrough = 0;

sizes.NumSampleTimes = 0;

sys = simsizes(sizes);

x0 =［0.1 0］;

str = ［］;

ts = ［］;

function sys=mdlDerivatives(t,x,u)

A=［0 1;0 -25］;

B=［0 133］';



ut=u(1);

dt=0.10*exp(-5*t);

 

dx=A*x+B*(ut+dt);

 

sys(1)=dx(1);

sys(2)=dx(2);

function sys=mdlOutputs(t,x,u)

sys(1)=x(1);

sys(2)=x(2);

(4) 控制器S函数： chap3_4ctrl.m

function ［sys,x0,str,ts］ = spacemodel(t,x,u,flag)

switch flag,

case 0,

［sys,x0,str,ts］=mdlInitializeSizes;

case 3,

sys=mdlOutputs(t,x,u);

case {2,4,9}

sys=［］;

otherwise

error(［'Unhandled flag = ',num2str(flag)］);

end

function ［sys,x0,str,ts］=mdlInitializeSizes

sizes = simsizes;

sizes.NumContStates= 0;

sizes.NumDiscStates  = 0;

sizes.NumOutputs     = 1;

sizes.NumInputs      = 2;

sizes.DirFeedthrough = 1;

sizes.NumSampleTimes = 1;

sys = simsizes(sizes);

x0  = ［］;

str = ［］;

ts  = ［0 0］;

function sys=mdlOutputs(t,x,u)

x1=u(1);

x2=u(2);

X=［x1 x2］';

 

xd=sin(t);dxd=cos(t);

ddxd=-sin(t);

 

e=x1-xd;de=x2-dxd;

X=［e de］';

 

%K =［-0.0099    0.1809］;

%K =［   -0.0305    0.1277］;

K =［  -0.0168    0.0667］;

 

tol=K*X;

ut=tol+1/133*ddxd+25/133*dxd;

 

sys(1)=ut;

(5) 作图程序： chap3_4plot.m

close all;

 

figure(1);

subplot(211);

plot(t,sin(t),'r',t,x(:,1),'b','linewidth',2);

xlabel('time(s)');ylabel('Angle tracking');

subplot(212);

plot(t,cos(t),'r',t,x(:,2),'b','linewidth',2);

xlabel('time(s)');ylabel('Angle speed tracking');

 

figure(2);

plot(t,ut(:,1),'r','linewidth',2);

xlabel('time(s)');ylabel('ut');

参考文献


［1］Grimm G, Hatfield J,Postlethwaite I, et al. Antiwindup for stable linear systems with input saturation: An LMIbased synthesis［J］. IEEE Transactions on Automatic Control, 2003,48(9): 15091525.

［2］Hu T S, Teel A, Zaccarian L. Stability and performance for saturated systems via quadratic and nonquadratic lyapunov functions［J］. IEEE Transactions on Automatic Control, 2006,51(3): 17701786.

［3］Wu H N,Zhu H Y,Wang J W. H∞ Fuzzy control for a class of nonlinear coupled ODEPDE systems with input constraint［J］. IEEE Transactions on Fuzzy Systems, 2015, 23(3): 393604.

［4］Fang H, Lin Z, Hu T. Analysis of linear systems in the presence of actuator saturation and L2disturbances［J］. Automatica, 2004, 40(7): 12291238.