T-S模糊模型与状态反馈控制及Matlab仿真 （附代码）

(一) 仿射非线性系统建模

以overhead crane system为例，

忽略摩擦，我们可以得到如下的状态空间方程，
x ˙ = ( x 2 m l x 4 2    sin ⁡    x 3 + m g    sin ⁡    x 3 cos ⁡    x 3    M + m    sin ⁡ 2    x 3 x 4 − m l x 4 2    sin ⁡    x 3    cos ⁡    x 3 + ( M + m ) g    sin ⁡    x 3    l ( M + m    sin ⁡ 2    x 3 ) ) + u M + m    sin ⁡ 2    x 3 ( 0 1 0 − cos ⁡    x 3 l ) \dot x=\begin{pmatrix}x_2\\\frac{mlx_4^2\;\sin\;x_3+mg\;\sin\;x_3\cos\;x_3\;}{M+m\;\sin^2\;x_3}\\x_4\\-\frac{mlx_4^2\;\sin\;x_3\;\cos\;x_3+(M+m)g\;\sin\;x_3\;}{l(M+m\;\sin^2\;x_3)}\end{pmatrix}+\frac u{M+m\;\sin^2\;x_3}\begin{pmatrix}0\\1\\0\\-\frac{\cos\;x_3}l\end{pmatrix} x˙=⎝⎜⎜⎜⎛​x2​M+msin2x3​mlx42​sinx3​+mgsinx3​cosx3​​x4​−l(M+msin2x3​)mlx42​sinx3​cosx3​+(M+m)gsinx3​​​⎠⎟⎟⎟⎞​+M+msin2x3​u​⎝⎜⎜⎛​010−lcosx3​​​⎠⎟⎟⎞​

(二) 计算T-S模糊模型子系统

取工作点。一般地，我们会选取平衡点。
x 0 = ( 0 0 0 0 ) x_0=\begin{pmatrix}0\\0\\0\\0\end{pmatrix} x0​=⎝⎜⎜⎛​0000​⎠⎟⎟⎞​, x 1 = ( 5 0 π / 6 π / 20 ) x_1=\begin{pmatrix}5\\0\\\pi/6\\\pi/20\end{pmatrix} x1​=⎝⎜⎜⎛​50π/6π/20​⎠⎟⎟⎞​, x 2 = ( 5 − 1 − π / 6 0 ) x_2=\begin{pmatrix}5\\-1\\-\pi/6\\0\end{pmatrix} x2​=⎝⎜⎜⎛​5−1−π/60​⎠⎟⎟⎞​
A 0 = J ( x 0 ) A_0=J(x_0) A0​=J(x0​), B 0 = g ( x 0 ) B_0=g(x_0) B0​=g(x0​)
A i = J ( x i ) + 1 x i T x i [ f ( x i ) − J ( x i ) x i ] x i T A_i=J(x_i)+\frac {1}{{x_i}^Tx_i}[f(x_i)-J(x_i)x_i]{x_i}^T Ai​=J(xi​)+xi​Txi​1​[f(xi​)−J(xi​)xi​]xi​T, B i = g ( x i ) B_i=g(x_i) Bi​=g(xi​), i = 1 , 2 i=1,2 i=1,2

(三) 建立推理，验证开环特性

if x is about x 0 x_0 x0​, then x ˙ = A 0 x + B 0 u \dot x=A_0x+B_0u x˙=A0​x+B0​u
if x is about x 1 x_1 x1​, then x ˙ = A 1 x + B 1 u \dot x=A_1x+B_1u x˙=A1​x+B1​u
if x is about x 2 x_2 x2​, then x ˙ = A 2 x + B 2 u \dot x=A_2x+B_2u x˙=A2​x+B2​u
取 h 0 ( x 3 ) = { sin ⁡    x 3 x 3 , x 3 ≠ 0 1 , x 3 = 0 h_0(x_3)=\left\{\begin{array}{l}\frac{\sin\;x_3}{x_3},x_3\neq0\\1,x_3=0\end{array}\right. h0​(x3​)={x3​sinx3​​,x3​​=01,x3​=0​， h 1 ( x 3 ) = h 2 ( x 3 ) = 1 − h 0 ( x 3 ) 2 h_1(x_3)=h_2(x_3)=\frac{1-h_0(x_3)}{2} h1​(x3​)=h2​(x3​)=21−h0​(x3​)​

(四) 极点配置，验证闭环特性

选取极点 p 0 = ( − 1.1 − 1.5 − 1.8 − 2.1 ) p_0=\begin{pmatrix}-1.1 -1.5 -1.8 -2.1\end{pmatrix} p0​=(−1.1−1.5−1.8−2.1​)， p 1 = ( − 1.3 − 1.6 − 1.9 − 2.2 ) p_1=\begin{pmatrix}-1.3 -1.6 -1.9 -2.2\end{pmatrix} p1​=(−1.3−1.6−1.9−2.2​)， p 2 = ( − 1.5 − 1.6 − 1.8 − 1.9 ) p_2=\begin{pmatrix}-1.5 -1.6 -1.8 -1.9\end{pmatrix} p2​=(−1.5−1.6−1.8−1.9​)

(五) 使用LMI验证稳定性

∃ P = ( 0.4016 0.1427 0.5521 0.5098 0.1427 0.1903 − 0.0489 0.7069 0.5521 − 0.0489 5.7520 − 0.1442 0.5098 0.7069 − 0.1442 2.8079 ) \exists P=\begin{pmatrix}0.4016&0.1427&0.5521&0.5098\\0.1427&0.1903&-0.0489&0.7069\\0.5521&-0.0489&5.7520&-0.1442\\0.5098&0.7069&-0.1442&2.8079\end{pmatrix} ∃P=⎝⎜⎜⎛​0.40160.14270.55210.5098​0.14270.1903−0.04890.7069​0.5521−0.04895.7520−0.1442​0.50980.7069−0.14422.8079​⎠⎟⎟⎞​
{ ( A i − B i K i ) T P + P ( A i − B i K i ) < 0 ,    ( i = 1 , 2 , 3 ) ( A i − B i K j + A j − B j K i ) T P + P ( A i − B i K j + A j − B j K i ) < 0 , ( i < j ≤ 3 ) \left\{\begin{array}{l}{(A_i-B_iK_i)}^TP+P(A_i-B_iK_i)<0,\;(i=1,2,3)\\{(A_i-B_iK_j+A_j-B_jK_i)}^TP+P(A_i-B_iK_j+A_j-B_jK_i)<0,(i<j\leq3)\end{array}\right. {(Ai​−Bi​Ki​)TP+P(Ai​−Bi​Ki​)<0,(i=1,2,3)(Ai​−Bi​Kj​+Aj​−Bj​Ki​)TP+P(Ai​−Bi​Kj​+Aj​−Bj​Ki​)<0,(i<j≤3)​

clear all;close all;
M=100;m=100;l=5;g=9.82;u=100;
syms x1 x2 x3 x4;
x00=[0 0 0 0]';x01=[5 0 pi/6 pi/20]';x02=[5 -1 -pi/6 0]';
x=[x1;x2;x3;x4];
f=[x2;(m*l*x4^2*sin(x3)+m*g*sin(x3)*cos(x3))/(M+m*sin(x3)*sin(x3));x4;-(m*l*x4^2*sin(x3)*cos(x3)+(M+m)*g*sin(x3))/l/(M+m*sin(x3)*sin(x3))];
gx=1/(M+m*sin(x3)*sin(x3))*[0;1;0;-cos(x3)/l];
J=jacobian(f,x);
Ai=J+(f-J*x)*x'/(x'*x);
A0=double(subs(J,x,x00));B0=double(subs(gx,x,x00));
A1=double(subs(Ai,x,x01));B1=double(subs(gx,x,x01));
A2=double(subs(Ai,x,x02));B2=double(subs(gx,x,x02));

a0=[-1.1 -1.5 -1.8 -2.1];
a1=[-1.3 -1.6 -1.9 -2.2];
a2=[-1.5 -1.6 -1.8 -1.9];

K0=place(A0,B0,a0);K1=place(A1,B1,a1);K2=place(A2,B2,a2);

setlmis([]);
P=lmivar(1,[4,1]);
lmiterm([1 1 1 P],1,A0-B0*K0,'s');
lmiterm([2 1 1 P],1,A1-B1*K1,'s');
lmiterm([3 1 1 P],1,A2-B2*K2,'s');
lmiterm([4 1 1 P],1,A0-B0*K1+A1-B1*K0,'s');
lmiterm([5 1 1 P],1,A1-B1*K2+A2-B2*K1,'s');
lmiterm([6 1 1 P],1,A2-B2*K0+A0-B0*K2,'s');
lmiterm([-7 1 1 P],1,1);
lmisys=getlmis;
[tmin,xfeas]=feasp(lmisys);
tmin
pmat=dec2mat(lmisys,xfeas,P)

t m i n = − 0.0163 t_{min}=-0.0163 tmin​=−0.0163
P = ( 4.2687 3.1372 6.7828 14.3290 3.1372 6.5107 − 10.3561 29.7178 6.7828 − 10.3561 151.4136 − 37.1458 14.3290 29.7178 − 37.1458 145.7880 ) P=\begin{pmatrix}4.2687&3.1372&6.7828&14.3290\\3.1372&6.5107&-10.3561&29.7178\\6.7828&-10.3561&151.4136&-37.1458\\14.3290&29.7178&-37.1458&145.7880\end{pmatrix} P=⎝⎜⎜⎛​4.26873.13726.782814.3290​3.13726.5107−10.356129.7178​6.7828−10.3561151.4136−37.1458​14.329029.7178−37.1458145.7880​⎠⎟⎟⎞​

下载链接
https://download.csdn.net/download/qq_29710939/84993530

trolly_sim.m

clear all;close all;
t0=0;tf=30;step=0.1;
x0(1:4)=[10 1 pi/10 0.1]';
global M m l u g;
global A0 A1 A2 B0 B1 B2 x3 K0 K1 K2;
M=100;m=50;l=4;g=9.81;u=100;
syms x1 x2 x3 x4;
x00=[0 0 0 0]';x01=[5 0 pi/6 pi/20]';x02=[5 -1 -pi/6 0]';
x=[x1;x2;x3;x4];
f=[x2;(m*l*x4^2*sin(x3)+m*g*sin(x3)*cos(x3))/(M+m*sin(x3)*sin(x3));x4;-(m*l*x4^2*sin(x3)*cos(x3)+(M+m)*g*sin(x3))/l/(M+m*sin(x3)*sin(x3))];
gx=1/(M+m*sin(x3)*sin(x3))*[0;1;0;-cos(x3)/l];
J=jacobian(f,x);
Ai=J+(f-J*x)*x'/(x'*x);
A0=double(subs(J,x,x00));B0=double(subs(gx,x,x00));
A1=double(subs(Ai,x,x01));B1=double(subs(gx,x,x01));
A2=double(subs(Ai,x,x02));B2=double(subs(gx,x,x02));

[t,x_fuzzy]=ode45('trolly_fuzzy',[t0:step:tf], x0);
[t,x_affine]=ode45('trolly_affine',[t0:step:tf], x0);
figure(1);
subplot(2,2,1);plot(t, x_fuzzy(:,1),'b-.',t,x_affine(:,1),'b:');xlabel('时间(秒)');ylabel('x1');
subplot(2,2,2);plot(t, x_fuzzy(:,2),'b-.',t,x_affine(:,2),'b:');xlabel('时间(秒)');ylabel('x2');
subplot(2,2,3);plot(t, x_fuzzy(:,3),'b-.',t,x_affine(:,3),'b:');xlabel('时间(秒)');ylabel('x3');
subplot(2,2,4);plot(t, x_fuzzy(:,4),'b-.',t,x_affine(:,4),'b:');xlabel('时间(秒)');ylabel('x4');

a0=[-1.1 -1.5 -1.8 -2.1];
a1=[-1.3 -1.6 -1.9 -2.2];
a2=[-1.5 -1.6 -1.8 -1.9];
K0=place(A0,B0,a0);K1=place(A1,B1,a1);K2=place(A2,B2,a2);

[t,x_fuzzy_ctrl]=ode45('trolly_fuzzy_ctrl',[t0:step:tf], x0);
[t,x_affine_ctrl]=ode45('trolly_affine_ctrl',[t0:step:tf], x0);
figure(2);
subplot(2,2,1);plot(t, x_fuzzy_ctrl(:,1),'b-.',t,x_affine_ctrl(:,1),'b:');xlabel('时间(秒)');ylabel('x1');
subplot(2,2,2);plot(t, x_fuzzy_ctrl(:,2),'b-.',t,x_affine_ctrl(:,2),'b:');xlabel('时间(秒)');ylabel('x2');
subplot(2,2,3);plot(t, x_fuzzy_ctrl(:,3),'b-.',t,x_affine_ctrl(:,3),'b:');xlabel('时间(秒)');ylabel('x3');
subplot(2,2,4);plot(t, x_fuzzy_ctrl(:,4),'b-.',t,x_affine_ctrl(:,4),'b:');xlabel('时间(秒)');ylabel('x4');

trolly_affine.m

function xdot = trolly_affine(t, x);
global M m l u g;
f=[x(2);(m*l*x(4)^2*sin(x(3))+m*g*sin(x(3))*cos(x(3)))/(M+m*sin(x(3))*sin(x(3)));x(4);-(m*l*x(4)^2*sin(x(3))*cos(x(3))+(M+m)*g*sin(x(3)))/l/(M+m*sin(x(3))*sin(x(3)))];
gx=1/(M+m*sin(x(3))*sin(x(3)))*[0;1;0;-cos(x(3))/l];
xdot=f+gx*u;

trolly_fuzzy.m

function xdot = trolly_fuzzy(t, x);
global A0 A1 A2 B0 B1 B2 u M m;
if(x(3)==0)
    h0=1;
else
    h0=sin(x(3))/x(3);
end
h1=(1-h0)/2;h2=h1;
xdot=h0*(A0*x+B0*u)+h1*(A1*x+B1*u)+h2*(A2*x+B2*u); 

trolly_affine_ctrl.m

function xdot = trolly_affine_ctrl(t, x);
global M m l u g K0 K1 K2;
f=[x(2);(m*l*x(4)^2*sin(x(3))+m*g*sin(x(3))*cos(x(3)))/(M+m*sin(x(3))*sin(x(3)));x(4);-(m*l*x(4)^2*sin(x(3))*cos(x(3))+(M+m)*g*sin(x(3)))/l/(M+m*sin(x(3))*sin(x(3)))];
gx=1/(M+m*sin(x(3))*sin(x(3)))*[0;1;0;-cos(x(3))/l];
h0=M*(sin(x(3)))^2/(M+m*sin(x(3))*sin(x(3)));
h1=M*(cos(x(3)))^2/(M+m*sin(x(3))*sin(x(3)));
h2=m*(sin(x(3)))^2/(M+m*sin(x(3))*sin(x(3)));
xdot=f+gx*(u-h0*K0*x-h1*K1*x-h2*K2*x);

trolly_fuzzy_ctrl.m

function xdot = trolly_fuzzy_ctrl(t, x);
global A0 A1 A2 B0 B1 B2 u M m K0 K1 K2;
if(x(3)==0)
    h0=1;
else
    h0=sin(x(3))/x(3);
end
h1=(1-h0)/2;h2=h1;
u0=u-h0*K0*x-h1*K1*x-h2*K2*x;
xdot=h0*(A0*x+B0*u0)+h1*(A1*x+B1*u0)+h2*(A2*x+B2*u0);