clear; clc; t_final = 26; tilde_u=[0.5;0.5]; x0=[-5;5]; %原始状态初始值 x=x0; %闭环初始状态 x_open=x0; %开环初始状态 % A1 = [0.9481 0;0 0.9481]; % A2 = [0.8279 0.032;0.049 0.29]; % B1 = [0.5 0.6;0.0562 0.0569]; % B2=B1; A1 = [1 0.1 ; -0.13 -0.1]; A2 = [1.1 0.15 ; -0.1 -0.2]; B1 = 13*[1.924 0.5124;0.1862 1.692]; B2=B1;% %系数矩阵B的维度 [n,m] = size(B1); C1 = [1.6 0;0 1.5]; C2 = C1; %对称正定权重矩阵 Q1 = [3 0;0 3]; Q2 = [3 0;0 3]; Q3 = [3 0;0 3]; R = 0.01*[0.05 0;0 0.05]; % R = 0.1*[0.1 0;0 1]; %设定初始值 y_tilde01 = [7;3]; %k-1时刻静态输出段的测量值，即K时刻经过FR协议的静态段输出值 y_tilde02 = [3;4]; %经过高速率FR协议后的动态输出段值 yinitial_static = [6;2]; yinitial_dynamic = [8;5];%测量输出初始值 y1=yinitial_static; y2=yinitial_dynamic; u_bar = 650; %输入约束 x_bar = 950; %状态约束 T=(C1)'*(inv(C1*(C1)')); t = 1:t_final; %RR协议 I = [1 0;0 1]; gamma = [1 0;0 0]; % gammab = I-gamma; %也就是文章中的T xk_open=zeros(2,t_final); xk = zeros(2,t_final); uk = zeros(2,t_final); Jk = zeros(1,t_final); y_tilde01_ks=zeros(2,t_final); y_tilde02_k=zeros(2,t_final); yk=zeros(2,t_final); modek = zeros(1,t_final); modek_static = zeros(1,t_final); gammak=zeros(1,t_final); gammak1=zeros(1,t_final); % V=[1 0;0 0]; V1=[1 0;0 0]; V2=[0 0;0 1]; %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%设定完初始值开始循环 %/ for k=t xk_open(:,k)=x_open; xk(:,k) = x; for ks=1:1:3 if mod(ks,2)==1 V=V1; end if mod(ks,2)==0 V=V2; end%此处for循环表达的应该是零阶保持器节点的选择（高速率情况下） tic % 表示开始的时间 setlmis([ ]); [Q_1p_un,nQ_1p_un,sQ_1p_un] = lmivar(1,[2,1]); %lmivar用来描述出现在线性矩阵不等式中的矩阵变量 [Q_2p_un,nQ_2p_un,sQ_2p_un] = lmivar(1,[2,1]);%可能是P [Q_3p_un,nQ_3p_un,sQ_3p_un] = lmivar(1,[2,1]); [Q_p_un,nQ_p_un,sQ_p_un] = lmivar(3,[sQ_1p_un, zeros(2,4);zeros(2),sQ_2p_un,zeros(2);zeros(2,4),sQ_3p_un]); [Q1c_co_un,n_Q1c_co_un,P1c_co_un]=lmivar(1,[2 1]);%Qr/Q_c [Q2c_co_un,n_Q2c_co_un,P2c_co_un]=lmivar(1,[2 1]); [Q3c_co_un,n_Q3c_co_un,P3c_co_un]=lmivar(1,[2 1]); [Q1z_co_un,n_Q1z_co_un,P1z_co_un]=lmivar(1,[2 1]);%Qt/Q_z [Q2z_co_un,n_Q2z_co_un,P2z_co_un]=lmivar(1,[2 1]); [Q3z_co_un,n_Q3z_co_un,P3z_co_un]=lmivar(1,[2 1]); [Q1cdot_co_un,n_Q1cdot_co_un,Q1c_dot_co_un]=lmivar(1,[2 1]);%Qr_head/dot [Q2cdot_co_un,n_Q2cdot_co_un,Q2c_dot_co_un]=lmivar(1,[2 1]); [Q3cdot_co_un,n_Q3cdot_co_un,Q3c_dot_co_un]=lmivar(1,[2 1]); [S11_un,nW1_un,sW1_un] = lmivar(1,[1,1]); [S22_un,nW2_un,sW2_un] = lmivar(1,[1,1]); [S_un,nW_un,sW_un] = lmivar(3,[sW1_un,0;0,sW2_un]);%[1 0;0 1]; tau_un = lmivar(1,[1,1]); %gamma标量/tau %所求的矩阵变量是正定的 Wlmi01=newlmi; lmiterm([-Wlmi01 1 1 Q1c_co_un],1,1); Wlmi02=newlmi; lmiterm([-Wlmi02 1 1 Q2c_co_un],1,1); Wlmi03=newlmi; lmiterm([-Wlmi03 1 1 Q3c_co_un],1,1); Wlmi04=newlmi; lmiterm([-Wlmi04 1 1 Q1z_co_un],1,1); Wlmi05=newlmi; lmiterm([-Wlmi05 1 1 Q2z_co_un],1,1); Wlmi06=newlmi; lmiterm([-Wlmi06 1 1 Q3z_co_un],1,1); Wlmi07=newlmi; lmiterm([-Wlmi07 1 1 Q1cdot_co_un],1,1); Wlmi08=newlmi; lmiterm([-Wlmi08 1 1 Q2cdot_co_un],1,1); Wlmi09=newlmi; lmiterm([-Wlmi09 1 1 Q3cdot_co_un],1,1); %第一个约束条件式17(无控制） lmi_017=newlmi; lmiterm([-lmi_017 1 1 S_un],T,1,'s'); %T代表T^l lmiterm([-lmi_017 1 1 Q1c_co_un],-1,1);%S11,S22,S都用S来表示 lmiterm([-lmi_017 2 2 S_un],1,1,'s'); lmiterm([-lmi_017 2 2 Q2c_co_un],-1,1,'s'); lmiterm([-lmi_017 3 3 S_un],1,1,'s'); lmiterm([-lmi_017 3 3 Q3c_co_un],-1,1);%第一个对角块表示的是Q_c^rightrowl lmiterm([-lmi_017 4 1 S_un],sqrtm(Q1)*T,1); lmiterm([-lmi_017 4 4 tau_un],1,1);%标量ri（单位矩阵）/tau lmiterm([-lmi_017 5 2 S_un],sqrtm(Q2),1); lmiterm([-lmi_017 5 5 tau_un],1,1); lmiterm([-lmi_017 6 3 S_un],sqrtm(Q3),1);% lmiterm([-lmi_017 6 6 tau_un],1,1);%第二个对角块表示的是phi_c^tilde lmiterm([-lmi_017 7 1 S_un],A1*T,1); lmiterm([-lmi_017 7 7 Q1z_co_un],1,1); lmiterm([-lmi_017 8 1 S_un],1,1);% lmiterm([-lmi_017 8 8 Q2z_co_un],1,1); lmiterm([-lmi_017 9 1 S_un],gamma,1);%C1代表的就是文章中C^[2]bar l，V不知道怎么变 lmiterm([-lmi_017 9 3 S_un],gammab,1);% lmiterm([-lmi_017 9 9 Q3z_co_un],1,1);%第三个对角块是psai^rightrow l，无控制，因此没有U代表的Y_tilde % lmi_032=newlmi; lmiterm([-lmi_032 1 1 0],eye(1)); lmiterm([-lmi_032 2 1 0],x_open);% lmiterm([-lmi_032 2 2 S_un],1,1); % lmi_033=newlmi; lmiterm([-lmi_033 1 1 0],2); lmiterm([-lmi_033 2 1 0],y_tilde01); lmiterm([-lmi_033 2 2 Q2c_co_un],1,1); lmiterm([-lmi_033 3 1 0],y_tilde02);%y_tilde0表示的是经过RR网络后的输出/FR（两者都是t_k-1时刻，且有关系表达式） lmiterm([-lmi_033 3 3 Q3c_co_un],1,1);%本文表述的是静态段和动态段两段的输出，两者表达没有联系/y^[1]（t_k-1) % lmi_035_0=newlmi; lmiterm([-lmi_035_0 1 1 Q1c_co_un],1,1); lmiterm([-lmi_035_0 1 1 Q1cdot_co_un],-1,1); lmi_035_1=newlmi; lmiterm([-lmi_035_1 1 1 Q1c_co_un],1,1); lmiterm([-lmi_035_1 1 1 S_un],-1,1); lmi_035_2=newlmi; lmiterm([-lmi_035_2 1 1 Q2c_co_un],1,1); lmiterm([-lmi_035_2 1 1 Q2cdot_co_un],-1,1); % lmi_0343=newlmi; lmiterm([-lmi_0343 1 1 Q3c_co_un],1,1); lmiterm([-lmi_0343 1 1 Q3cdot_co_un],-1,1); % lmi_44=newlmi; % lmiterm([-lmi_44 1 1 0],eye(2)); % lmiterm([-lmi_44 2 1 X_Yr_co_un'],1,1); % lmiterm([-lmi_44 2 2 S_un],u_bar*u_bar,1); % lmi_048=newlmi; lmiterm([-lmi_048 1 1 0],eye(2)); lmiterm([-lmi_048 2 1 Q1c_co_un],1,1); lmiterm([-lmi_048 2 2 Q1c_co_un],x_bar*x_bar,1); plwsys=getlmis;%为什么这个线性矩阵不等式的变量名是一样的？？？？？？？？ options=[1e-5,0,0,0,0]; c=[1 zeros(1,38)]; [gamma_coopt_un,xopt]=mincx(plwsys,c,options); P1r_coopt_un=dec2mat(plwsys,xopt,Q1c_co_un); P2r_coopt_un=dec2mat(plwsys,xopt,Q2c_co_un); P3r_coopt_un=dec2mat(plwsys,xopt,Q3c_co_un); P1t_coopt_un=dec2mat(plwsys,xopt,Q1z_co_un); P2t_coopt_un=dec2mat(plwsys,xopt,Q2z_co_un); P3t_coopt_un=dec2mat(plwsys,xopt,Q3z_co_un); Q1r_coopt_un=dec2mat(plwsys,xopt,Q1cdot_co_un); Q2r_coopt_un=dec2mat(plwsys,xopt,Q2cdot_co_un); Q3r_coopt_un=dec2mat(plwsys,xopt,Q3cdot_co_un); Q_popt_un = dec2mat(plwsys,xopt,Q_p_un); S_popt_un = dec2mat(plwsys,xopt,S_un); %%%%%%有控制 setlmis([ ]); [Q_1p,nQ_1p,sQ_1p] = lmivar(1,[2,1]); %lmivar用来描述出现在线性矩阵不等式中的矩阵变量 [Q_2p,nQ_2p,sQ_2p] = lmivar(1,[2,1]); [Q_3p,nQ_3p,sQ_3p] = lmivar(1,[2,1]); [Q_p,nQ_p,sQ_p] = lmivar(3,[sQ_1p, zeros(2,4);zeros(2),sQ_2p,zeros(2);zeros(2,4),sQ_3p]); [Q_1c_co,n_Q1c_co,P1c_co]=lmivar(1,[2 1]);%Qr（tilde） [Q_2c_co,n_Q2c_co,P2c_co]=lmivar(1,[2 1]); [Q_3c_co,n_Q3c_co,P3c_co]=lmivar(1,[2 1]); [Q_1z_co,n_Q1z_co,P1z_co]=lmivar(1,[2 1]);%Qt（tilde） [Q_2z_co,n_Q2z_co,P2z_co]=lmivar(1,[2 1]); [Q_3z_co,n_Q3z_co,P3z_co]=lmivar(1,[2 1]); [Q1cdot_co,n_Q1cdot_co,Q1c_dot_co]=lmivar(1,[2 1]);%Qr_hat [Q2cdot_co,n_Q2cdot_co,Q2c_dot_co]=lmivar(1,[2 1]); [Q3cdot_co,n_Q3cdot_co,Q3c_dot_co]=lmivar(1,[2 1]); [X_Y2c_co,N_Y2c_co,Y2c_co]=lmivar(2,[2 2]);%Y2r [X_N2c_co,N_N2c_co,N2c_co]=lmivar(2,[2 2]);% [S11,nW1,sW1] = lmivar(1,[1,1]); % [S22,nW2,sW2] = lmivar(1,[1,1]); [S,nW,sW] = lmivar(3,[sW1,0;0,sW2]); tau = lmivar(1,[1,1]); %标量 %所求的矩阵变量是正定的 Wlmi1=newlmi; lmiterm([-Wlmi1 1 1 Q_1c_co],1,1); Wlmi2=newlmi; lmiterm([-Wlmi2 1 1 Q_2c_co],1,1); Wlmi3=newlmi; lmiterm([-Wlmi3 1 1 Q_3c_co],1,1); Wlmi4=newlmi; lmiterm([-Wlmi4 1 1 Q_1z_co],1,1); Wlmi5=newlmi; lmiterm([-Wlmi5 1 1 Q_2z_co],1,1); Wlmi6=newlmi; lmiterm([-Wlmi6 1 1 Q_3z_co],1,1); Wlmi7=newlmi; lmiterm([-Wlmi7 1 1 Q1cdot_co],1,1); Wlmi8=newlmi; lmiterm([-Wlmi8 1 1 Q2cdot_co],1,1); Wlmi9=newlmi; lmiterm([-Wlmi9 1 1 Q3cdot_co],1,1); %第一个约束条件式17(有控制) lmi_17=newlmi; lmiterm([-lmi_17 1 1 S],T,1,'s'); lmiterm([-lmi_17 1 1 Q_1c_co],-1,1); lmiterm([-lmi_17 2 2 S],1,1,'s'); lmiterm([-lmi_17 2 2 Q_2c_co],-1,1); lmiterm([-lmi_17 3 3 S],1,1,'s'); lmiterm([-lmi_17 3 3 Q_3c_co],-1,1); lmiterm([-lmi_17 4 1 S],sqrtm(Q1)*T,1); lmiterm([-lmi_17 4 4 tau],1,1); lmiterm([-lmi_17 5 2 S],sqrtm(Q2),1); lmiterm([-lmi_17 5 5 tau],1,1); lmiterm([-lmi_17 6 3 S],sqrtm(Q3),1);% lmiterm([-lmi_17 6 6 tau],1,1); lmiterm([-lmi_17 7 1 X_Y2c_co],sqrtm(R)*gamma,1);% lmiterm([-lmi_17 7 2 X_Y2c_co],sqrtm(R),1);% lmiterm([-lmi_17 7 3 X_Y2c_co],sqrtm(R)*gammab,1);% lmiterm([-lmi_17 7 7 tau],1,1); lmiterm([-lmi_17 8 1 S],A1*T,1); lmiterm([-lmi_17 8 1 X_Y2c_co],B1*gamma,1);%V lmiterm([-lmi_17 8 2 X_Y2c_co],B1,1);%v1 lmiterm([-lmi_17 8 3 X_Y2c_co],B1*gammab,1);%v3 lmiterm([-lmi_17 8 8 Q_1z_co],1,1); lmiterm([-lmi_17 9 1 S],1,1);% lmiterm([-lmi_17 9 9 Q_2z_co],1,1); lmiterm([-lmi_17 10 1 S],gamma,1);% lmiterm([-lmi_17 10 3 S],gammab,1);% lmiterm([-lmi_17 10 10 Q_3z_co],1,1); % lmi_32=newlmi; lmiterm([-lmi_32 1 1 0],eye(1)); lmiterm([-lmi_32 2 1 0],x); lmiterm([-lmi_32 2 2 S],1,1); % lmi_33=newlmi; lmiterm([-lmi_33 1 1 0],2); lmiterm([-lmi_33 2 1 0],y_tilde01);%y[1] lmiterm([-lmi_33 2 2 Q_2c_co],1,1); lmiterm([-lmi_33 3 1 0],y_tilde02);%y[2]bar lmiterm([-lmi_33 3 3 Q_3c_co],1,1); % lmi_35_0=newlmi; lmiterm([-lmi_35_0 1 1 Q_1c_co],1,1); lmiterm([-lmi_35_0 1 1 Q1cdot_co],-1,1); lmi_35_1=newlmi; lmiterm([-lmi_35_1 1 1 Q_1c_co],1,1); lmiterm([-lmi_35_1 1 1 S],-1,1); lmi_35_2=newlmi; lmiterm([-lmi_35_2 1 1 Q_2c_co],1,1); lmiterm([-lmi_35_2 1 1 Q2cdot_co],-1,1); lmi_35_3=newlmi; lmiterm([-lmi_35_3 1 1 Q_3c_co],1,1); lmiterm([-lmi_35_3 1 1 Q3cdot_co],-1,1); % lmi_47=newlmi; lmiterm([-lmi_47 1 1 0],eye(2)); lmiterm([-lmi_47 2 1 X_Y2c_co'],1,1); %Yr对应Y2r lmiterm([-lmi_47 2 2 S11],u_bar*u_bar,1); %这个s对应的是s11 % lmi_48=newlmi; lmiterm([-lmi_48 1 1 0],eye(2)); lmiterm([-lmi_48 2 1 Q_1c_co],1,1); %048表示的是没控制的 lmiterm([-lmi_48 2 2 Q_1c_co],x_bar*x_bar,1); plwsys=getlmis; options=[1e-5,0,0,0,0]; c=[1 zeros(1,46)]; [gamma_coopt,xopt]=mincx(plwsys,c,options); P1c_coopt=dec2mat(plwsys,xopt,Q_1c_co); P2c_coopt=dec2mat(plwsys,xopt,Q_2c_co); P3c_coopt=dec2mat(plwsys,xopt,Q_3c_co); P1z_coopt=dec2mat(plwsys,xopt,Q_1z_co); P2z_coopt=dec2mat(plwsys,xopt,Q_2z_co); P3z_coopt=dec2mat(plwsys,xopt,Q_3z_co); Q1c_coopt=dec2mat(plwsys,xopt,Q1cdot_co); Q2c_coopt=dec2mat(plwsys,xopt,Q2cdot_co); Q3c_coopt=dec2mat(plwsys,xopt,Q3cdot_co); Y2c_coopt=dec2mat(plwsys,xopt,X_Y2c_co); Q_popt = dec2mat(plwsys,xopt,Q_p);%未知 S_popt = dec2mat(plwsys,xopt,S); N2C_popt = dec2mat(plwsys,xopt,X_N2c_co); %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%矩阵变量的最优解，线性矩阵不等式的表达与赋值已完成 Fc = Y2c_coopt/S_popt;%Q_r,Q_t(tilde),Q_r(hat),F_r,Y_2r,S11, lambda1 = rand(1); lambda2 = 1-lambda1; A = lambda1*A1+lambda2*A2; B = lambda1*B1+lambda2*B2; C_1 = lambda1*C1+lambda2*C2; C_2 = lambda1*C1+lambda2*C2; y_bar =y_tilde01+ gamma*C_2*x + gammab*y_tilde02;%公式七 （kelt...函数，ZOH的part）（y_tilde表示的是经过高速率和RR网络协议后的值） u = Fc*y_bar;%公式九（输出反馈，控制律，控制增益F_c） %开环系统 x_open_temp=A*x_open; %闭环系统 x_temp = A*x+ B*u;%公式1 y2 = C_2*x; y1 = C_1*x; %更新状态值 x_open=x_open_temp;%temp表示更新值 x=x_temp; %y[1] y_tilde01_ks(:,ks)=y_tilde01; y_tilde01_ks_temp=C_1*x; y_tilde01=y_tilde01_ks_temp; %ybar[2] y_tilde02_k(:,k)=y_tilde02; y_tilde02_k_temp=gamma*C_2*x + gammab*y_tilde02; y_tilde02=y_tilde02_k_temp; %%存一下值 yk(:,ks)=y1; uk(:,k)=u; gammak(:,k)=gamma_coopt; % gammak1(:,k)=gamma_coopt_un; if (gamma(1,1)==1) mode = 1; elseif(gamma(1,1)~=1) mode = 2; end modek(:,k) = mode; end%这块是生成figure3的 theta1 = (y2(1,1)-y_tilde02(1,1))^2; theta2 = (y2(2,1)-y_tilde02(2,1))^2; if (theta1 >= theta2) gamma = [1 0;0 0]; elseif (theta1 < theta2) gamma = [0 0; 0 1];% end gammab = eye(2) - gamma;%更新传输通道 end %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %开环闭环系统对比 figure(1); subplot(2,1,1) plot(t-1,xk_open(1,:),'b-.',t-1,xk(1,:),'r--'); legend('the open-loop system','the closed-loop system'); xlabel('Time(k)'), ylabel('x1'); subplot(2,1,2) plot(t-1,xk_open(2,:),'b-.',t-1,xk(2,:),'r--'); legend('the open-loop system','the closed-loop system'); xlabel('Time(k)'), ylabel('x2'); %时间间隔d不同时的状态对比 % figure(2); % subplot(2,1,1) % plot(t-1,xk(1,:),'c-o'); % legend('d=3'); % xlabel('Time(k)'),ylabel('State value x_1' ,'FontSize',10,'Rotation',90); % subplot(2,1,2) % plot(t-1,xk(2,:),'m-x'); % legend('d=3'); % xlabel('Time(k)'),ylabel('State value x_2' ,'FontSize',10,'Rotation',90); %信号转换 % figure (3); % stairs(modek,'b','LineWidth',1.25); % axis([0 t_final 0 3]) % set(gca,'ytick', [0 1 2 3]) % legend({'Switching signal $$\vartheta(k)$$'},'Interpreter','latex'); % xlabel({'Time step $$(k)$$'},'Interpreter','latex'); % ylabel({'Switching signal'},'Interpreter','latex'); % figure(4); % subplot(2,1,1) % plot(t-1,gammak(:,t),'r-'); % xlabel('Time(k)'),ylabel('\gamma'); % subplot(2,1,2) % plot(t-1,gammak1(:,t),'r-'); % xlabel('Time(k)'),ylabel('\gammak1'); % % figure(5); % plot (t-1,yk(:,t),'b-'); % xlabel('Time(k)'),ylabel('\gamma'); % figure(6); subplot(2,1,1) plot (t-1,xk(1,:),'m-',t-1,xk(2,:),'b--'); legend('State value x_1','State value x_2'); xlabel('Time(k)'),ylabel('x'); subplot(2,1,2) plot (t-1,uk(1,:),'m-',t-1,uk(2,:),'b--'); legend('control inpuut u1','control inpuut u2'); xlabel('Time(k)'),ylabel('u');