Find Q and R matriks to get eigenvalue positive
2 views (last 30 days)
Show older comments
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)
In this code, i need to find the s1,s2,s3,s4,s7,s12 to get my matrix S.
After i get the S matrix, I find the eigenvalue of matrix 4x4
I have problem to get the true Q and R matrix, if I can get the true Q and R matrix I can get the eigenvalue positive.
0 Comments
Answers (0)
See Also
Categories
Find more on Linear Algebra in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!