Teljes Matlab script kiegészítő függvényekkel.
File: d2018_01_13_Kalman_decomp_LPV_generate_matrices.m Directory: projects/3_outsel/2018_01_10_LPV_inversion Author: Peter Polcz (ppolcz@gmail.com)
Created on 2018. January 13.
Generate matrices for:
LPV decomposition (version 2)
File: d2018_01_13_Kalman_decomp_LPV.m Directory: projects/3_outsel/2018_01_10_LPV_inversion Author: Peter Polcz (ppolcz@gmail.com)
Created on 2018. January 13.
Given the following LPV system:
$$ \begin{aligned} &\dot x = A(\rho) x + B(\rho) u,~~~ \rho \in \mathcal P \\ &y = C x \\ &\begin{aligned} \text{where: } & A(\rho) = A_0 + A_1 \rho \in \mathbb{R}^{n\times n} \\ & B(\rho) = B_0 + B_1 \rho \in \mathbb{R}^{n\times r} \\ & C \in \mathbb{R}^{m\times n} \\ & D = 0_{m\times r} \end{aligned} \end{aligned} $$
First, generate matrices $A_0$, $B_0$, $C$, such that $(A_0, B_0, C)$ be uncontrollable and unobservable.
try c = evalin('caller','persist'); catch; c = []; end
persist = pcz_persist(mfilename('fullpath'), c); clear c;
persist.backup();
%clear persist
│ - Persistence for `d2018_01_13_Kalman_decomp_LPV_generate_matrices` initialized [run ID: 48, 318]Persistence for `2018.01.14. Sunday, 00:37:56` │ - Script `d2018_01_13_Kalman_decomp_LPV_generate_matrices` backuped
tol = 1e-10;
a = 2; % Nr. of contr. and obs.
b = 2; % Nr. of contr. and unobs.
c = 1; % Nr. of uncontr. and obs.
d = 2; % Nr. of uncontr. and unobs.
n = a+b+c+d; % Nr. of states
r = 2; % Nr. of inputs
m = 2; % Nr. of outputs
A_ = [
randn(a,a) zeros(a,b) randn(a,c) zeros(a,d)
randn(b,a) randn(b,b) randn(b,c) randn(b,d)
zeros(c,a) zeros(c,b) randn(c,c) zeros(c,d)
zeros(d,a) zeros(d,b) randn(d,c) randn(d,d)
];
B_ = [
randn(a+b,round(r/2)) * randn(round(r/2),r)
% randn(a,r)
% randn(b,r)
zeros(c,r)
zeros(d,r)
];
C_ = [
randn(m,a) zeros(m,b) randn(m,c) zeros(m,d)
];
D0 = zeros(m,r);
System $(\bar A, \bar B, \bar C, D)$ is in Kalman decomposed form. The transfer function is reducible exactly to an $a$th order system.
fprintf '\nSize of (A_,B_,C_,D0): \n'
sys = ss(A_,B_,C_,D0);
size(minreal(sys,tol,false))
Size of (A_,B_,C_,D0): State-space model with 2 outputs, 2 inputs, and 2 states.
Transform the system with a random transformation matrix $T$.
Random permutation matrix
T = pcz_permat(randperm(n));
Random "nice" transformation matrix
while 1
T = round(randn(n)) .* (rand(n) > 0.8);
if abs(det(T)) < 0.9 || abs(det(T)) > 20 || sum(sum(T ~= 0)) > 2*n
continue
end
Ti = inv(T);
if sum(sum(Ti ~= 0)) > 2*n
continue
end
break
end
% pcz_display(T)
A0 = T \ A_ * T;
B0 = T \ B_;
C0 = C_ * T;
fprintf '\nSize of (A0,B0,C0,D0): \n'
sys = minreal(ss(A0,B0,C0,D0),tol,false);
size(sys)
Size of (A0,B0,C0,D0): State-space model with 2 outputs, 2 inputs, and 2 states.
Now, generate $A(\rho) = A_0 + A_1 \rho$ with a sparse matrix $A_1$. The same with $B$.
A1 = randn(n).*(rand(n) > 0.95);
B1 = randn(n,r).*(rand(n,m) > 0.94);
C1 = zeros(m,n);
D1 = zeros(m,r);
pcz_display(A1,B1);
A1 [7x7] =
0 0 0 0 0.4885 0 0
0 0 -0.6896 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 0 0 0
0 0 0 0 0 0 0
B1 [7x2] =
0 0
0 0
0 0
0 0
0 0
0 0
0 0
Check joint controllability and observability for a few values of $\rho$.
A_fh = @(rho) A0 + A1*rho;
B_fh = @(rho) B0 + B1*rho;
C_fh = @(rho) C0 + C1*rho;
D_fh = @(rho) D0 + D1*rho;
H = @(rho) minreal(ss(A_fh(rho), B_fh(rho), C_fh(rho), D_fh(rho)),tol,false);
fprintf '\nSize of (A(ρ),B(ρ),C(ρ),D(ρ)): \n'
for rho = -1:0.5:1
% minreal(ss(A_fh(rho), B_fh(rho), C_fh(rho), D_fh(rho)))
size(H(rho))
end
Size of (A(ρ),B(ρ),C(ρ),D(ρ)): State-space model with 2 outputs, 2 inputs, and 4 states. State-space model with 2 outputs, 2 inputs, and 4 states. State-space model with 2 outputs, 2 inputs, and 2 states. State-space model with 2 outputs, 2 inputs, and 4 states. State-space model with 2 outputs, 2 inputs, and 4 states.
A = {A0, A1};
B = {B0, B1};
C = {C0, C1};
D = {D0, D1};
AA = [A{:}];
BB = [B{:}];
tol = 1e-10;
Im_B = pcz_vecalg_orth(BB,tol);
Ker_C = pcz_vecalg_ints(pcz_vecalg_null(C0,tol), pcz_vecalg_null(C1,tol));
[R,V] = CSA(AA,Im_B,Ker_C);
% Check the strong invertibility condition, the followings has to be empty
Strong_invertibility = pcz_vecalg_ints(V,Im_B,tol)'
return
Strong_invertibility = 0×7 empty double matrix
Export values
pcz_display(a)
pcz_display(b)
pcz_display(c)
pcz_display(d)
pcz_display(n)
pcz_display(r)
pcz_display(m)
pcz_num2str_multiline(A_,'format','%8.4g')
pcz_num2str_multiline(B_,'format','%8.4g')
pcz_num2str_multiline(C_,'format','%8.4g')
pcz_num2str_multiline(T,'format','%3.g')
pcz_num2str_multiline(A0,'format','%8.4g')
pcz_num2str_multiline(A1,'format','%8.4g')
pcz_num2str_multiline(B0,'format','%8.4g')
pcz_num2str_multiline(B1,'format','%8.4g')
pcz_num2str_multiline(C0,'format','%8.4g')
pcz_num2str_multiline(C1,'format','%1g')
pcz_num2str_multiline(D0,'format','%1g')
pcz_num2str_multiline(D1,'format','%1g')
cd([proot 'projects/3_outsel/2018_01_10_LPV_inversion']);
save(['LPV_' persist.date], '-regexp', '^[a,b,c,d,n,m,r]$', ...
'^[A,B,C]_$', '^[A,B,C,D]_fh$', ...
'^T$','^[A,B,C,D]0$','^[A,B,C,D]1$')