Find Q and R matriks to get eigenvalue positive

Question

Ivan Dwi Putra 2020-6-8

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/544166-find-q-and-r-matriks-to-get-eigenvalue-positive

My code

clear; close; clc;
syms a1_head a2_head b hstar 
%Parameter Massa
m1 = 8095; % massa train set 1 dalam kg
m2 = 8500; % massa train set 2 dalam kg
g  = 10;
c_0_1 = 0.01176;
c_1_1 = 0.00077616;
c_2_1 = 4.48 ;
c_0_2 = 0.01176 ;
c_1_2 = 0.00077616;
c_2_2 = 4.48;
v_0 = 300;
hstar = 120;
a2_head
a_1   = -1./m1.*(c_1_1 + 2.*c_2_1.*v_0);  
a_2   = -1./m2.*(c_1_2 + 2.*c_2_2.*v_0); 
a_1_head = 1-(a_1.*hstar);
a_2_head = 1-(a_2.*hstar);
b = 1;
% Model data
A = sym(zeros(4,4));
A(1,2) = a_1_head;
A(3,2) = (a_2_head) - 1;  A(3,4) = a_2_head;
display(A);
B = sym(zeros(4,2));
B(1,1) = -b*hstar;
B(2,1) = b;
B(3,2) = -b*hstar   ;
B(4,1) = -b; B(4,2) = b;
display(B);
% Q and R matrices for ARE
Q = sym(zeros(4,4)); Q(1,:) = [10000 100 100 100]; Q(2,:) = [100 10000 100 100]; Q(3,:) = [100 100 10000 100]; Q(4,:) = [100 100 100 10000]; display(Q);
R = sym(zeros(2,2)); R(1,:) = [2500 90]; R(2,:) = [90 5000]; display(R);
% S matrix Value
s1_1 =  -(124919^(1/2)*50000000^(1/2))/60000;
s1_2 =   (124919^(1/2)*50000000^(1/2))/60000;
s2_1 = -5.2761e+03;
s2_2 =  5.2761e+03;
s3_1 = 13.4620;
s3_2 = 16.1414;
s4_1 = -232.2135;
s4_2 = 322.5401;
s7_1 = -1.5834e+05;
s7_2 =  1.5522e+05;
s12_1 = 19.1007;
s12_2 = 4.8186;
S_matrix_1 = [s1_1 s2_1 s3_1  s4_1;
              s2_1 s1_1 s7_1  s2_1;
              s3_1 s7_1 s1_1 s12_1;
              s4_1 s2_1 s12_1 s1_1
            ];
          
S_matrix_2 = [s1_2 s2_2 s3_2  s4_2;
              s2_2 s1_2 s7_2  s2_2;
              s3_2 s7_2 s1_2 s12_2;
              s4_2 s2_2 s12_2 s1_2
            ];
        
eig_1 = eig(S_matrix_1);
eig_2 = eig(S_matrix_2);
K1 = inv(R)*transpose(B)*S_matrix_1;
display(K1);
K2 = inv(R)*transpose(B)*S_matrix_2;
display(K2);
% Matrix S to find
svar = sym('s',[1 16]);
S = [svar(1:4); svar(5:8); svar(9:12); svar(13:16)];
S(2,1) = svar(2);
S(2,2) = svar(1);
S(2,4) = svar(2);
S(3,1) = svar(3);
S(3,2) = svar(7);
S(3,3) = svar(1);
S(4,1) = svar(4);
S(4,2) = svar(2);
S(4,3) = svar(12);
S(4,4) = svar(1);
display(S);
% LHS of ARE: A'*S + S*A' - S*B*Rinv*B'*S
left_ARE = transpose(A)*S + S*A - S*B*inv(R)*transpose(B)*S; 
display(left_ARE);
% RHS of ARE: A'*S + S*A' - S*B*Rinv*B'*S
right_ARE = -Q;
display(right_ARE);
% Find S in ARE
C = left_ARE;
display(C);
D = right_ARE;
display(D);
% Element Matriks S
C1 = C(1,1);
display(C1);
C2 = C(1,2);
display(C2);
C3 = C(1,3);
display(C3);
C4 = C(1,4);
display(C4);
C5 = C(2,1);
display(C5);
C6 = C(2,2);
display(C6);
C7 = C(2,3);
display(C7);
C8 = C(2,4);
display(C8);
C9 = C(3,1);
display(C9);
C10 = C(3,2);
display(C10);
C11 = C(3,3);
display(C11);
C12 = C(3,4);
display(C12);
C13 = C(4,1);
display(C13);
C14 = C(4,2);
display(C14);
C15 = C(4,3);
display(C15);
C16 = C(4,4);
display(C16);
D1 = D(1,1);
display(D1);
D2 = D(1,2);
display(D2);
D3 = D(1,3);
display(D3);
D4 = D(1,4);
display(D4);
D5 = D(2,1);
display(D5);
D6 = D(2,2);
display(D6);
D7 = D(2,3);
display(D7);
D8 = D(2,4);
display(D8);
D9 = D(3,1);
display(D9);
D10 = D(3,2);
display(D10);
D11 = D(3,3);
display(D11);
D12 = D(3,4);
display(D12);
D13 = D(4,1);
display(D13);
D14 = D(4,2);
display(D14);
D15 = D(4,3);
display(D15);
D16 = D(4,4);
display(D16);
E1 = C1 == D1;
display(E1);
E2 = C2 == D2;
display(E2);
E3 = C3 == D3;
display(E3)
E4 = C4 == D4;
display(E4);
E5 = C5 == D5;
display(E5);
E6 = C6 == D6;
display(E6);
E7 = C7 == D7;
display(E7);
E8 = C8 == D8;
display(E8);
E9 = C9 == D9;
display(E9);
E10 = C10 == D10;
display(E10);
E11 = C11 == D11;
display(E11);
E12 = C12 == D12;
display(E12);
E13 = C13 == D13;
display(E13);
E14 = C14 == D14;
display(E14);
E15 = C15 == D15;
display(E15);
E16 = C16 == D16;
display(E16);
syms s1 s2 s3 s4 s5 s6 s7 s8 s9 s10 s11 s12 s13 s14 s15 s16
[Sol_s1_1, Sol_s2_1, Sol_s3_1, Sol_s4_1] = solve(E1,s1,s2,s3,s4)
[Sol_s2_2, Sol_s3_2, Sol_s4_2, Sol_s7_1] = solve(E2,s2,s3,s4,s7,'Real',true)
[Sol_s2_3, Sol_s3_3, Sol_s4_3, Sol_s7_2, sol_s12_1] = solve(E3,s2,s3,s4,s7,s12,'Real',true) 
[Sol_s3_4, Sol_s4_4,sol_s12_2] = solve(E4,s3,s4,s12,'Real',true) 
[Sol_s4_5, Sol_s7_3] = solve(E5,s4,s7,'Real',true)
[Sol_s7_4] = solve(E6,s7,'Real',true)
  
[Sol_s7_5, Sol_s12_3] = solve(E7,s7,s12,'Real',true)
[Sol_s12_4] = solve(E8,s12,'Real',true)
% 
% [Sol_s2_9, Sol_s3_7, Sol_s4_7, Sol_s7_7, Sol_s12_5] = solve(E9,s2,s3,s4,s7,s12,'Real',true)
% 
% [Sol_s3_8, Sol_s7_8, Sol_s12_6] = solve(E10,s3,s7,s12,'Real',true)
% 
% [Sol_s3_9, Sol_s7_9, Sol_s12_7] = solve(E11,s3,s7,s12,'Real',true)
% 
% [Sol_s4_8, Sol_s7_10, Sol_s12_8] = solve(E12,s4,s7,s12,'Real',true)
% 
% [Sol_s4_9, Sol_s12_9] = solve(E13,s4,s12,'Real',true)
% 
% [Sol_s7_11, Sol_s12_10] = solve(E14,s7,s12,'Real',true)
% 
% [Sol_s12_11] = solve(E15,s12,'Real',true)
% 
% [Sol_s12_12] = solve(E16,s12,'Real',true)
% syms s1 s2 s3 s4 s5 s6 s7 s8 s9 s10 s11 s12 s13 s14 s15 s16
% [Sol_s1_1, Sol_s2_1, Sol_s3_1, Sol_s4_1] = solve(E1,s1,s2,s3,s4)
% 
% [Sol_s2_2, Sol_s3_2, Sol_s4_2, Sol_s7_1] = solve(E2,s2,s3,s4,s7,'Real',true)
%  
% [Sol_s3_3, Sol_s4_3, Sol_s7_2, sol_s12_1] = solve(E3,s3,s4,s7,s12,'Real',true) 
%  
% [Sol_s4_4, Sol_s12_2] = solve(E4,s4,s12,'Real',true) 
% 
% [Sol_s7_3] = solve(E5,s7,'Real',true)
% % 
% % [Sol_s8_4] = solve(E6,s8,'Real',true)
% %   
% [Sol_s12_3] = solve(E7,s12)
%  
% % [Sol_s8_6, Sol_s12_4, Sol_s16] = solve(E8,s8,s12,s16,'Real',true)
% %  
% % [Sol_s11_5, Sol_s12_5] = solve(E9,s11,s12,'Real',true)
% %  
% % [Sol_s12_6] = solve(E10,s12,'Real',true)
% %  
% % [Sol_s12_7] = solve(E11,s12,'Real',true)
% % 
% % [Sol_s16_3] = solve(E12,s16,'Real',true)
% % % % 
% % % [Sol_s8_9, Sol_s12_9, Sol_s16_4] = solve(E13,s8,s12,s16,'Real',true)
% % % 
% % % [Sol_s2_22, Sol_s4_10, Sol_s6_7, Sol_s7_11, Sol_s8_10, Sol_s12_10, Sol_s16_5] = solve(E14,s2,s4,s6,s7,s8,s12,s16,'Real',true)
% % % 
% % % [Sol_s3_12, Sol_s4_11, Sol_s7_12, Sol_s8_11, Sol_s11_9, Sol_s12_11, Sol_s16_6] = solve(E15,s3,s4,s7,s8,s11,s12,s16,'Real',true)
% % % 
% % % [Sol_s4_12, Sol_s8_12, Sol_s12_12, Sol_s16_7] = solve(E16,s4,s8,s12,s16,'Real',true)