5x5 matrix (P) with all 25 values unknown. I have two known matrices of size 5x5 each (A1 and A2) and the relation P*A1*inv(P)=A2. How do I obtain P?

10 次查看(过去 30 天)
A2=[1.1581 0.0279 -0.0739 -0.0169 -0.0100
-0.9172 0.7637 0.5688 0.1299 0.0795
12.8599 2.5037 -5.4500 -1.4719 -0.8824
-39.7814 -7.8095 20.0738 5.5809 2.7482
-27.9052 -5.5369 14.1918 3.2387 2.9446];
A1=[ -1.0093 0.0645 -0.0795 -0.0172 -0.5886
4.1505 0.0040 -0.6850 -0.0824 -2.9194
1.1250 0.7191 -2.4753 -0.0098 -0.3402
-1.1250 1.1521 -0.2040 -1.2556 -0.4159
2.5296 0.0455 -0.1234 -0.0646 -3.7925];
syms P a1 b1 c1 d1 e1 a2 b2 c2 d2 e2 a3 b3 c3 d3 e3 a4 b4 c4 d4 e4 a5 b5 c5 d5 e5
P1=[a1 b1 c1 d1 e1
a2 b2 c2 d2 e2
a3 b3 c3 d3 e3
a4 b4 c4 d4 e4
a5 b5 c5 d5 e5];
V=vpa(P1*A1*inv(P1),6);
solve(V==A2)
%%This gives all 0s.
%%I also tried the below code but still getting zeroes.
global A1 A2 B1 B2
guess=[10;0.1; 0.1; 0.1; 0.1; 0.1; 10; 0.1; 0.1; 0.1; 0.1; 0.1; 10; 0.1; 0.1; 0.1; 0.1; 0.1; 1; 0.1; 0.1; 0.1; 0.1; 0.1; 1];
%guess=[1 -1 -2 1 2 1 1 -3 2 4 -2 1 1 3 2 3 -3 1 1.5 2.5 2.8 3 3.5 3.2 1];
%guess=[1 0.1 0.1 0.1 0.1 0.1 1 0.1 0.1 0.1 0.1 0.1 1 0.1 0.1 0.1 0.1 0.1 1 0.1 0.1 0.1 0.1 0.1 1];
options = optimoptions(@lsqnonlin,'OptimalityTolerance', 1e-16, 'FunctionTolerance', 1e-16,'StepTolerance',1e-16);
result=lsqnonlin(@eqns,guess,[],[],options);
%result=fsolve(@eqns,guess,options);
P=[result(1) result(6) result(11) result(16) result(21)
result(2) result(7) result(12) result(17) result(22)
result(3) result(8) result(13) result(18) result(23)
result(4) result(9) result(14) result(19) result(24)
result(5) result(10) result(15) result(20) result(25)];
function fncs = eqns(z)
global A1 A2 B1 B2
a1=z(1);
a2=z(2);
a3=z(3);
a4=z(4);
a5=z(5);
b1=z(6);
b2=z(7);
b3=z(8);
b4=z(9);
b5=z(10);
c1=z(11);
c2=z(12);
c3=z(13);
c4=z(14);
c5=z(15);
d1=z(16);
d2=z(17);
d3=z(18);
d4=z(19);
d5=z(20);
e1=z(21);
e2=z(22);
e3=z(23);
e4=z(24);
e5=z(25);
P1=[a1 b1 c1 d1 e1
a2 b2 c2 d2 e2
a3 b3 c3 d3 e3
a4 b4 c4 d4 e4
a5 b5 c5 d5 e5];
expres=P1*A1*inv(P1)-A2;
fncs(1)= expres(1,1);
fncs(2)= expres(1,2);
fncs(3)= expres(1,3);
fncs(4)= expres(1,4);
fncs(5)= expres(1,5);
fncs(6)= expres(2,1);
fncs(7)= expres(2,2);
fncs(8)= expres(2,3);
fncs(9)= expres(2,4);
fncs(10)= expres(2,5);
fncs(11)= expres(3,1);
fncs(12)= expres(3,2);
fncs(13)= expres(3,3);
fncs(14)= expres(3,4);
fncs(15)= expres(3,5);
fncs(16)= expres(4,1);
fncs(17)= expres(4,2);
fncs(18)= expres(4,3);
fncs(19)= expres(4,4);
fncs(20)= expres(4,5);
fncs(21)= expres(5,1);
fncs(22)= expres(5,2);
fncs(23)= expres(5,3);
fncs(24)= expres(5,4);
fncs(25)= expres(5,5);
fncs(26)= expres(4,5);
end
Even if we assume that A1 and A2 do not satisfy this relation A2=P*A1*inv(P), how do we find P for some other A1 A2 which do satisfy this relation. Please help!
  2 个评论
Richa Dubey
Richa Dubey 2022-1-22
Lets take A1 and A2 both as Identity matrices of order 5. Then P=eye(5) shall be the exact value of P which satisfies A2=P*A1*inv(P). How do we obtain this P by solving this equation in matlab code?
Matt J
Matt J 2022-1-22
Then P=eye(5) shall be the exact value of P which satisfies A2=P*A1*inv(P).
That will be one solution, but clearly there will be infinitely many alternative solutions as well. Any invertible P would solve the equation.

请先登录,再进行评论。

采纳的回答

Matt J
Matt J 2022-1-22
编辑:Matt J 2022-1-22
A necessary but not sufficient condition for a solution is A2*P-P*A1=0, which is equivalent to K*P(:)=0 where K=kron(I,A2)-kron(A1.',I). The necessary condition can therefore be solved (if a solution exists) with null(). Example:
A1=[1.1581 0.0279 -0.0739 -0.0169 -0.0100
-0.9172 0.7637 0.5688 0.1299 0.0795
12.8599 2.5037 -5.4500 -1.4719 -0.8824
-39.7814 -7.8095 20.0738 5.5809 2.7482
-27.9052 -5.5369 14.1918 3.2387 2.9446];
P0=rand(5);
A2=P0*A1/P0;
d=length(A1);
I=eye(d);
K=kron(I,A2)-kron(A1.',I);
P=reshape(null(K),d,d,[]);
errors=nan(size(P,3),1);
for j=1:numel(errors)
Pj=P(:,:,j);
errors(j)=max(abs(A2-Pj*A1/Pj),[],'all');
end
P=P(:,:,errors<1e-8)
P =
P(:,:,1) = -0.1484 -0.0809 0.3049 0.1126 -0.0748 -0.0825 0.4494 -0.5312 -0.1467 -0.2100 -0.1645 0.1687 -0.0116 0.0257 -0.1725 -0.1118 0.2680 -0.2137 -0.0442 -0.1774 -0.1682 0.1182 -0.0494 -0.0055 -0.0969 P(:,:,2) = 0.1499 0.1345 0.1784 0.0311 0.2248 0.3233 0.1632 -0.2085 0.2133 0.1304 0.3195 0.1930 -0.0885 0.0778 0.2600 0.2801 0.2077 -0.1852 0.1193 0.1836 0.4238 0.1047 -0.0101 0.0908 0.1241 P(:,:,3) = 0.2255 0.3731 -0.4906 -0.0203 -0.1859 -0.2177 0.0766 0.1463 0.1861 -0.1393 0.0645 0.3138 -0.2742 0.0879 -0.2131 -0.0570 0.2065 -0.0975 0.1219 -0.1739 0.0441 0.1757 -0.1544 0.0525 -0.1212 P(:,:,4) = 0.1178 -0.0384 0.2374 0.1743 -0.2159 0.2590 0.0815 0.2378 0.1002 -0.2551 0.2093 0.0079 0.3701 0.2168 -0.3097 0.2308 0.0235 0.3178 0.1624 -0.2681 0.0283 0.0140 0.1750 0.0982 -0.1668 P(:,:,5) = 0.0787 -0.1273 -0.2835 -0.0928 -0.1639 0.4456 -0.0677 -0.2022 -0.3277 -0.0984 0.2488 -0.1367 -0.2216 -0.1928 -0.1933 0.3397 -0.1347 -0.1554 -0.2293 -0.1354 -0.0409 -0.0693 -0.1722 -0.1530 -0.0908
This does produce several solutions, though it is not clear to me if solving the necessary conditions will always find one.
  4 个评论
Matt J
Matt J 2022-1-22
编辑:Matt J 2022-1-22
This method isnt working when null(K) in the above code returns an empty matrix.
If null(K) is empty, then the method has succeeded. It has successfully determined that no solutions exists. For a solution to exist, it must lie in the nullspace of K and therefore null(K) cannot be empty.
Could you please suggest a soln which is 'necessary and sufficient' for satisfying A2=P*A1*inv(P)?
See my comment above about fminunc.
Richa Dubey
Richa Dubey 2022-1-23
Thank you so much!
The soln you gave works well but only for those A1 and A2 which satisfy this relation exactly.
Although I was wondering how to deal with those A1 and A2 which don't satisfy A2=P*A1*inv(P) exactly, as in the value of A2 (LHS) and P*A1*inv(P) (RHS) differ by a small error of order say 10^(-3).
How do we calculate P for such A1 and A2 so that this error is minimized. I tried dealing with this using 'lsqnonlin' for a 3rd orded system though but wasn't satisfied with the results.
guess=[1 0.5 0.5 0.5 1 0.5 0.5 0.5 1];
options = optimoptions(@lsqnonlin,'OptimalityTolerance', 1e-100,'Algorithm','levenberg-marquardt', 'FunctionTolerance', 1e-100,'StepTolerance',1e-100,'MaxIterations',1e+100,'MaxFunctionEvaluations',1e+50,'FiniteDifferenceType','central','ScaleProblem','jacobian');
result=lsqnonlin(@eqns,guess,[],[],options);
%result=fsolve(@eqns,guess,options);
P=[result(1) result(4) result(7)
result(2) result(5) result(8)
result(3) result(6) result(9)]
function fncs = eqns(z)
global A1 A2 B1 B2
a1=z(1);
a2=z(2);
a3=z(3);
b1=z(4);
b2=z(5);
b3=z(6);
c1=z(7);
c2=z(8);
c3=z(9);
P1=[a1 b1 c1
a2 b2 c2
a3 b3 c3];
expres=P1*A1*inv(P1)-A2;
expres1=P1*B1-B2;
fncs(1)= expres(1,1);
fncs(2)= expres(1,2);
fncs(3)= expres(1,3);
fncs(4)= expres(2,1);
fncs(5)= expres(2,2);
fncs(6)= expres(2,3);
fncs(7)= expres(3,1);
fncs(8)= expres(3,2);
fncs(9)= expres(3,3);
fncs(10)= expres1(1);
fncs(11)= expres1(2);
fncs(12)= expres1(3);
end
I calculated error using this:
A22=P*A1*inv(P);
B22=P*B1;
eA=((A22-A2)./A2)*100
eB=((B22-B2)./B2)*100

请先登录,再进行评论。

更多回答(1 个)

Torsten
Torsten 2022-1-22
编辑:Torsten 2022-1-22
I think you could also get such P by solving the linear optimization problem
min: sum_ij (eij+ + eij-)
under the constraints
E+ - E- = (P+ - P-)*A1 - A2*(P+ - P-)
sum_ij (pij+ + pij-) = 1
E+,E-,P+,P- >= 0
The matrix P = P+ - P- should have the property P*A1 = A2*P if the optimal value of the objective is 0.
Otherwise, such P does not exist.
Note that it's not guaranteed that P is invertible (I guess). But it's different from the trivial solution P=0.
  1 个评论
Richa Dubey
Richa Dubey 2022-1-23
Thanks for the soln.
I am unable to implement it though. Could you please elaborate more on this?
Also, please look into my comment in the above ans if possible.

请先登录,再进行评论。

类别

Help CenterFile Exchange 中查找有关 Robust Control Toolbox 的更多信息

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by