一般线性多智能体系统的包容控制

多智能体包容控制containment control

连续时间下静态控制器的仿真如下：

李忠奎老师11年的一篇文章，题目如下

Distributed containment control of multi-agent systems with general linear dynamics in the presence of multiple leaders
智能体动力学模型

$$\begin{array}{l}{\dot{x}}_i=A{x}_i+B{u}_i,\\ y_i=C{x}_i,i=1,\cdot \cdot \cdot ,N,\end{array}$$

控制协议

$$\begin{array}{l}{\dot{v}}_i=A{v}_i+B{u}_i+cL\sum _{j\in F\cup R}{a}_{ij}\left[C\left(v_i-v_j\right)-\left(y_i-y_j\right)\right],\\ u_i=cK\sum _{j\in F\cup R}{a}_{ij}\left(v_i-v_j\right),i\in F\end{array}$$

%example 1:     continuous-time static containment controller
%code author:   hzwrabbit

%% initial state
clc;clear;close all;
global A B L1 L2 nf nl
nf=6;
nl=3;
A=[0 0    0       1      0      0;
   0 0    0       0      1      0;
   0 0    0       0      0      1;
   0 0 -0.2003 -0.2003   0      0;
   0 0  0.2003    0   -0.2003   0;
   0 0    0       0      0    -1.6129];
B=[ 0       0;
    0       0;
    0       0;
  0.9441 0.9441;
  0.9441 0.9441;
-28.7079 28.7079];

L1=[ 3  0  0 -1 -1 -1 ;
   -1  1  0  0  0  0 ;
   -1 -1  2  0  0  0 ;
   -1  0  0  2  0  0 ;
    0  0  0 -1  2  0;
    0  0  0  0 -1  2];
L2=[0  0  0;
    0  0  0;
    0  0  0;
    0  0 -1;
    0 -1  0;
   -1  0  0];

L=[ 3  0  0 -1 -1 -1  0  0  0;
   -1  1  0  0  0  0  0  0  0;
   -1 -1  2  0  0  0  0  0  0;
   -1  0  0  2  0  0  0  0 -1;
    0  0  0 -1  2  0  0 -1  0;
    0  0  0  0 -1  2 -1  0  0;
    0  0  0  0  0  0  0  0  0;
    0  0  0  0  0  0  0  0  0;
    0  0  0  0  0  0  0  0  0];
[V,D]=eig(L1);
c_th=1/min(diag(D));        % c_th=1.2176, c>=c_th,取c=2

%% LMI求解P
setlmis([])
P=lmivar(1,[6 1]);
lmiterm([1 1 1 P],1,A','s');
lmiterm([-1 1 1 0],2*B*B');
lmiterm([-2 1 1 P],1,1);
lmis1 = getlmis;
[tmin,xfeas] = feasp(lmis1);
%判断是否有可行解
if (tmin < 0)  
    disp('Feasible');
    P = dec2mat(lmis1,xfeas,P);
else
    P = nan;
end
global F
F=-B'*inv(P);

% F=[-0.0089, 0.0068, 0.0389,-0.0329, 0.0180, 0.0538;
%     0.0068,-0.0089,-0.0389, 0.0180,-0.0329,-0.0538];

%% calculate ODE funtion
tspan=[0,50];
x=2*(rand(54,1)-0.5);
[t,y]=ode45(@contain,tspan,x);

%% draw picture

figure(1)

for i=1:(nf+nl-1)
    if(i<6)
        plot(t,y(:,1+nf*i),'b','linestyle',':','linewidth',1.5);
    else
        plot(t,y(:,1+nf*i),'r','linewidth',2);
    end
grid on;
hold on;
end
xlabel('t');
ylabel('x_{i1}');

figure(2)
for i=1:(nf+nl-1)
    if(i<6)
        plot(t,y(:,2+nf*i),'b','linestyle',':','linewidth',1);
    else
        plot(t,y(:,2+nf*i),'r','linewidth',1.5);
    end
grid on;
hold on;
end
xlabel('t');
ylabel('x_{i2}');

figure(3)
for i=1:(nf+nl-1)
    if(i<6)
        plot(t,y(:,3+nf*i),'b','linestyle',':','linewidth',1);
    else
        plot(t,y(:,3+nf*i),'r','linewidth',1.5);
    end
grid on;
hold on;
end
xlabel('t');
ylabel('x_{i3}');

%% function 
function   ydot=contain(t,x)
c=2;  
global A B F L1 L2
a11=kron(eye(6),A)+c*kron(L1,B*F);
a12=c*kron(L2,B*F);
a21=zeros(18,36);
a22=kron(eye(3),A);
ydot=[a11 a12;
      a21 a22]*x;
end`

初值由rand函数生成，仿真如下

4.2离散时间下包容控制器

仿真时遇到的问题：计算文章中给的D1的特征值时与matlab不同(经过线上matlab的验算确认)

按文章所给特征值算的delta求取参数矩阵K和L与文章相同。

最后仿真结果和文章不同

clc;clear;close all;
%% initial state
global K L A B C D1 D2  nf nl K L
nf=6;
nl=3;
A=[1 1;
   0 1];
B=[0 1]';
C=[1 0];
D1=[ 0.4  0    0  0.1   0.3 0.2;
     0.5 0.5   0   0     0   0;
     0.3 0.2  0.5  0     0   0;
     0.5  0    0  0.5    0   0;
      0   0    0  0.4   0.4 0.2;
      0   0    0   0    0.3 0.7];
  [V,X]=eig(D1)

D2=[0  0  0;
    0  0  0;
    0  0  0;
    0  0  0;
    0  0  0;
    0  0  0];

D=[ 0.4  0    0  0.1   0.3 0.2   0  0   0;
     0.5 0.5   0   0     0   0   0  0   0;
     0.3 0.2  0.5  0     0   0   0  0   0;
     0.5  0    0  0.5    0   0   0  0   0;
      0   0    0  0.4   0.4 0.2  0  0   0;
      0   0    0   0    0.3 0.7  0  0   0
      0   0    0   0     0   0   1  0   0;
      0   0    0   0     0   0   0  1   0;    
      0   0    0   0     0   0   0  0   1];   
      
 

% c_th=max(diag(X))
Lamda_max=0.9;
delta=1 - Lamda_max*Lamda_max;
%% 描述修正离散代数riccati方程求K
Q=[1 0;
   0 1];
P1=care(A,B,Q);
re=A'*P1*A - delta*A'*P1*B*inv(B'*P1*B)*B'*P1*A + Q - P1;
% 迭代求解修正离散代数riccati方程
e=1;
while(e>1e-10)             %达到精度要求
    Pe=A'*P1*A - delta*A'*P1*B*inv(B'*P1*B)*B'*P1*A + Q ;
    e=norm(Pe-P1);
    P1=Pe;
end
P1;
K=inv(B'*P1*B)*B'*P1*A;

%% 描述修正离散代数riccati方程求L

P2=care(A,B,Q);
re=A*P2*A' - delta*A*P2*C'*inv(C*P2*C')*C*P2*A' + Q - P2;
% re=A'*P2*A - delta*A'*P2*C'*inv(C*P2*C')*C*P2*A + R - P2;
% 迭代求解修正离散代数riccati方程
e=1;
while(e>1e-10)             %达到精度要求
    Pe=A*P2*A' - delta*A*P2*C'*inv(C*P2*C')*C*P2*A' + Q;
    e=norm(Pe-P2);
    P2=Pe;
end
P2
L=-A*P2*C'*inv(C*P2*C')

%% calculate ODE funtion
tspan=[0,10];

x0=4*(rand(36,1)-0.5);

[t,y]=ode45(@contain,tspan,x0);

%% draw picture
figure(1)
% plot(t,y(:,1),t,y(:,5),t,y(:,9),t,y(:,13),t,y(:,17),'b','linestyle',':','linewidth',1.5);
% plot(t,y(:,25),t,y(:,29),t,y(:,33),'linewidth',2.5,'color','r');
for i=1:8
    if(i<6)
        plot(t,y(:,1+4*i),'b','linestyle',':','linewidth',1.5);
    else
        plot(t,y(:,1+4*i),'r','linewidth',1.8);
    end
grid on;
hold on;
end
xlabel('t');
ylabel('x_{i1}');
% figure(1)
% plot(t,y(:,1),t,y(:,7),t,y(:,13),t,y(:,19),t,y(:,25),t,y(:,31),'linewidth','0.8');
% % plot(t,y(:,37),t,y(:,43),t,y(:,49),'linewidth',3.5);
% plpt(t,y(:,37),t,y(:,43),t,y(:,49),'linewidth','3.5','color','r');

figure(2)
for i=1:8
    if(i<6)
        plot(t,y(:,2+4*i),'b','linestyle',':','linewidth',1);
    else
        plot(t,y(:,2+4*i),'r','linewidth',1.5);
    end
grid on;
hold on;
end
xlabel('t');
ylabel('x_{i2}');


%% function 
function   ydot=contain(t,x)
global K L A B  C D1 D2  nf nl K L
s11=A;
s12=zeros(length(A));
s21=zeros(length(A));
s22=A;
S=[s11 s12;
   s21 s22];
h11=zeros(length(A));
h12=B*K;
h21=-L*C;
h22=L*C+B*K;
H=[h11 h12;
   h21 h22];

a11=kron(eye(nf),S)+kron((eye(nf)-D1),H);
a12=-kron(D2,H);
a21=zeros(12,24);
a22=kron(eye(nl),S);
ydot=[a11 a12;
      a21 a22]*x;
end

仿真的example2中，算出D1的特征值和文中不一样，做出的仿真的也有问题，做出的朋友可以讨论讨论。本人QQ号：1911936192