Teljes Matlab script kiegészítő függvényekkel.
File: Kalman_decomposition.m Directory: demonstrations/oktatas/ccs/2017fall Author: Peter Polcz (ppolcz@gmail.com)
Created on 2018. January 13.
Inhereted from:
CCS 2017 fall. Homework 2. Solutions
File: d2017_10_26_hf2_mo.m Directory: demonstrations/oktatas/ccs/2017fall Author: Peter Polcz (ppolcz@gmail.com)
Created on 2017. October 04. Reviewed on 2017. October 30.
% Automatically generated stuff
global SCOPE_DEPTH VERBOSE
SCOPE_DEPTH = 0;
VERBOSE = 1;
a = 4; % 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,r)
randn(b,r)
zeros(c,r)
zeros(d,r)
];
C_ = [
randn(m,a) zeros(m,b) randn(m,c) zeros(m,d)
];
D = zeros(m,r);
System $(\bar A, \bar B, \bar C, D)$ is in Kalman decomposed form. The transfer function is reducible to a $a$rd order system.
sys = ss(A_,B_,C_,D)
sys =
A =
x1 x2 x3 x4 x5 x6
x1 2.18 -2.601 -1.168 0.3713 0 0
x2 0.08162 1.055 -2.237 0.8478 0 0
x3 0.6023 1.042 -0.007656 -0.5526 0 0
x4 0.3258 0.7171 -0.8829 -0.4567 0 0
x5 0.4452 0.1199 -0.2424 -0.4134 -1.39 2.35
x6 -0.04099 -0.4354 1.042 -1.472 0.869 -1.258
x7 0 0 0 0 0 0
x8 0 0 0 0 0 0
x9 0 0 0 0 0 0
x7 x8 x9
x1 0.3771 0 0
x2 0.1719 0 0
x3 -1.037 0 0
x4 -1.203 0 0
x5 -0.1934 -1.313 -0.3133
x6 0.9227 1.149 0.1707
x7 0.02033 0 0
x8 -1.577 0.324 0.4185
x9 1.445 -1.152 -0.5577
B =
u1 u2
x1 0.4645 -0.5897
x2 0.5284 -0.103
x3 -0.3536 -1.79
x4 0.5544 -0.2775
x5 -0.17 -0.01351
x6 0.4434 -0.4604
x7 0 0
x8 0 0
x9 0 0
C =
x1 x2 x3 x4 x5 x6 x7
y1 1.939 0.5564 1.077 -1.643 0 0 -0.181
y2 -0.6088 -2.241 -0.1898 -1.077 0 0 -0.413
x8 x9
y1 0 0
y2 0 0
D =
u1 u2
y1 0 0
y2 0 0
Continuous-time state-space model.
Transfer function
H = minreal(tf(sys))
H =
From input 1 to output...
-0.09687 s^3 + 1.399 s^2 - 10.95 s - 4.647
1: ---------------------------------------------
s^4 - 2.771 s^3 + 2.831 s^2 - 3.084 s - 5.319
-1.997 s^3 + 0.482 s^2 + 7.199 s + 12.67
2: ---------------------------------------------
s^4 - 2.771 s^3 + 2.831 s^2 - 3.084 s - 5.319
From input 2 to output...
-2.673 s^3 + 8.612 s^2 - 20.4 s - 20.13
1: ---------------------------------------------
s^4 - 2.771 s^3 + 2.831 s^2 - 3.084 s - 5.319
1.228 s^3 - 13.58 s^2 + 19.27 s + 24.21
2: ---------------------------------------------
s^4 - 2.771 s^3 + 2.831 s^2 - 3.084 s - 5.319
Continuous-time transfer function.
Transform the system with a random orthogonal transformation matrix $T$.
T = orth(rand(n));
A = T * A_ * T';
B = T * B_;
C = C_ * T';
sys = ss(A,B,C,D)
sys =
A =
x1 x2 x3 x4 x5 x6 x7
x1 0.6481 0.3399 0.2673 -0.1689 0.7096 0.2047 -0.9442
x2 1.472 -1.058 1.253 -0.2828 0.06571 0.9535 0.01389
x3 -1.11 2.086 0.7939 0.415 0.9133 0.07849 -1.042
x4 0.2161 -0.7022 1.778 0.878 0.2339 0.4193 0.5418
x5 0.3258 0.2156 -0.3072 1.648 0.8113 -0.7606 -0.3288
x6 -0.3634 0.1001 -0.5182 1.44 0.8654 -0.9298 0.8704
x7 0.5611 0.6293 -0.6845 0.01012 0.6761 -0.3199 -0.6573
x8 -0.05946 0.5881 1.097 0.3105 1.252 0.1918 -1.436
x9 1.48 -0.8469 -0.2679 0.02202 0.4993 0.4339 -0.2813
x8 x9
x1 1.367 0.1473
x2 1.171 -0.2404
x3 -0.6024 0.2386
x4 0.09242 -0.03678
x5 -0.3055 0.1739
x6 0.9128 0.4748
x7 1.179 0.211
x8 0.3095 -1.015
x9 0.7015 -0.8854
B =
u1 u2
x1 -0.624 0.07024
x2 0.1063 0.7054
x3 0.2961 0.6281
x4 0.21 0.5862
x5 0.02116 -0.108
x6 -0.5239 0.9498
x7 -0.4097 0.6003
x8 -0.344 -1.145
x9 0.2437 -0.1614
C =
x1 x2 x3 x4 x5 x6 x7
y1 0.3165 -0.9533 -1.394 -0.5288 -0.08648 -0.55 -1.576
y2 0.9369 0.7279 -1.228 -0.9609 -0.5572 1.057 1.183
x8 x9
y1 -0.3537 -1.346
y2 -0.2139 0.2145
D =
u1 u2
y1 0 0
y2 0 0
Continuous-time state-space model.
Controllability matrix
Cn = ctrb(A,B);
Transformation matrix, which generates the controllability staircase form.
[S1,~,~] = svd(Cn)
S1 =
Columns 1 through 7
-0.0344 0.3801 0.3272 0.2280 0.0698 0.5390 -0.2868
-0.6202 0.2357 0.1980 -0.2393 0.0155 -0.0890 -0.5261
0.2950 0.4295 -0.4772 -0.3318 -0.0880 -0.3618 -0.3844
-0.4784 0.3112 -0.3427 -0.2697 0.2311 0.1248 0.5737
0.2629 0.3668 -0.0962 -0.1104 -0.5778 0.4586 0.1417
0.1937 0.4314 0.2660 -0.0783 0.4344 -0.0199 0.2114
0.2279 0.3233 0.4762 0.0706 0.0093 -0.4843 0.2097
-0.1725 0.2991 -0.3651 0.8274 -0.0094 -0.1772 -0.0200
-0.3316 0.0623 0.2617 0.0026 -0.6411 -0.2801 0.2345
Columns 8 through 9
-0.5636 0.0001
0.3995 -0.1370
-0.3009 0.1170
-0.2129 -0.2019
0.4258 -0.1711
0.3293 0.5977
-0.0596 -0.5742
0.1746 0.0227
-0.2512 0.4584
The transformation.
A1 = S1' * A * S1
B1 = S1' * B
C1 = C * S1
A1 =
Columns 1 through 7
-2.3836 -1.6226 0.5723 -1.0918 -1.1167 -0.1776 0.3619
-1.1182 2.0772 -1.9653 -0.0905 -1.0177 1.1393 -0.0918
-0.0468 0.0820 0.7525 2.8086 -0.0779 0.7580 0.4491
0.0492 -0.0970 -0.8506 0.3153 0.0427 1.0901 -0.7856
0.0001 0.0000 -0.0003 0.0000 -0.4074 -1.6492 0.2009
-0.0000 0.0000 0.0002 -0.0001 -0.1051 -0.2308 0.6883
0.0000 0.0000 -0.0000 -0.0000 0.0000 -0.0000 1.2414
0.0000 -0.0000 -0.0000 -0.0000 0.0000 -0.0000 0.9952
0.0000 -0.0000 0.0000 -0.0000 0.0000 0.0000 0.6748
Columns 8 through 9
1.4716 -0.9392
0.5885 -0.1770
-0.5567 -0.1298
0.8485 0.2451
0.1743 -0.1406
-0.8109 -0.1991
-0.6782 -1.1084
-0.7094 -0.8543
0.1764 -0.7454
B1 =
-0.2683 0.0083
-0.4580 0.8570
-0.5436 0.5866
-0.5968 -1.4871
-0.4160 0.6909
-0.2155 -0.2901
0.0000 -0.0000
0.0000 0.0000
-0.0000 -0.0000
C1 =
Columns 1 through 7
0.4408 -1.8361 -0.3503 0.5506 0.6704 1.8680 -0.1243
-0.0929 0.2844 2.3984 0.5919 0.6191 -0.1071 -0.2836
Columns 8 through 9
0.1252 0.0402
0.2858 0.0918
Order of the controllable subsystem
ab = rank(Cn);
2.1. Observability staircase form of the controllable subsystem
On1 = obsv(A1(1:ab,1:ab), C1(:,1:ab));
[~,~,S21] = svd(On1);
2.2. Observability staircase form of the uncontrollable subsystem
On2 = obsv(A1(ab+1:end,ab+1:end), C1(:,ab+1:end));
[~,~,S22] = svd(On2);
Second transformation matrix to generate the observability staircase form for both controllable and uncontrollable subsystems.
S2 = blkdiag(S21, S22)
S2 =
Columns 1 through 7
-0.1520 -0.0154 -0.1615 0.0525 0.0207 0.9733 0
0.6509 0.1025 0.6771 -0.1870 -0.1417 0.2287 0
-0.4077 -0.5483 0.3817 -0.5462 0.2981 0.0141 0
-0.5960 0.4149 0.5743 0.3436 -0.1572 -0.0064 0
-0.1620 -0.0507 -0.1320 -0.4055 -0.8884 -0.0072 0
-0.0747 0.7168 -0.1500 -0.6176 0.2769 0.0022 0
0 0 0 0 0 0 0.6868
0 0 0 0 0 0 -0.6920
0 0 0 0 0 0 -0.2222
Columns 8 through 9
0 0
0 0
0 0
0 0
0 0
0 0
0.4034 -0.6046
0.6173 -0.3743
-0.6754 -0.7032
Transformation
A2 = round(S2' * A1 * S2,10)
B2 = round(S2' * B1,10)
C2 = round(C1 * S2,10)
A2 =
Columns 1 through 7
2.3794 0.3914 0.0509 0.0015 0 0 0.1258
1.3450 0.0269 -1.0693 0.0081 0 0 -0.1514
0.8970 3.0656 1.1008 -0.0040 0 0 -0.7720
0.5112 -0.1264 -1.4905 -0.7363 0 0 -1.4350
-1.0939 0.8399 0.2354 -0.9976 0.0945 0.0230 0.7058
0.4528 -0.4343 -0.7251 0.1151 1.5036 -2.7422 -0.6249
0 0 0 0 0 0 0.0203
0 0 0 0 0 0 1.1324
0 0 0 0 0 0 -1.8152
Columns 8 through 9
0 0
0 0
0 0
0 0
-0.3333 0.0998
1.7421 -0.1143
0 0
0.0714 0.2435
-1.3266 -0.3051
B2 =
0.4034 1.1135
-0.1257 -1.0940
-0.7298 -0.0989
0.4651 -1.0923
0.3010 -0.4067
-0.3673 0.2163
0 0
0 0
0 0
C2 =
Columns 1 through 7
-1.6956 1.5306 -1.5007 -0.6786 0 0 -0.1810
-1.2236 -1.1471 1.3973 -1.3496 0 0 -0.4130
Columns 8 through 9
0 0
0 0
Using minreal the transfer function which generates the Kalman decomposition can be retained easily, however, this does not work all the times. The uncontrollable subsystem is not always decomposed (since the uncontrollable subsystem is already redundant).
[H,U] = minreal(sys);
A1 = round(U*A/U,10)
B1 = round(U*B,10)
C1 = round(C/U,10)
5 states removed.
A1 =
Columns 1 through 7
2.4045 -0.7869 -1.3048 0.7232 0 0 0.5497
1.7864 1.6418 0.8533 1.7592 0 0 0.4637
-0.6565 -1.4494 -0.4373 0.6298 0 0 1.2209
0.0515 -0.9723 -0.3470 -0.8381 0 0 -0.8284
-0.8031 0.3060 -0.0494 0.9086 0.1067 -1.4807 -0.3325
0.6820 -1.1515 -0.3126 -0.6556 0 -2.7543 0.8821
0 0 0 0 0 0 0.0203
0 0 0 0 0 0 0.6476
0 0 0 0 0 0 2.0391
Columns 8 through 9
0 0
0 0
0 0
0 0
-0.5109 -0.0945
-1.6069 -0.5631
0 0
-0.2116 -0.2196
1.3506 -0.0222
B1 =
0.4155 0.1923
-0.8576 -0.3107
0.1212 -1.8310
-0.0672 -0.3915
0.0948 -0.2587
-0.4653 0.3811
0 0
0 0
0 0
C1 =
Columns 1 through 7
-1.3675 -0.5879 0.9465 2.1951 0 0 0.1810
-1.8037 1.3804 -0.9356 -0.7434 0 0 0.4130
Columns 8 through 9
0 0
0 0
dim_X_C = rank(Cn);
X_C = S1(:,1:dim_X_C)
Alternatively
X_C = orth(Cn);
X_C =
-0.0344 0.3801 0.3272 0.2280 0.0698 0.5390
-0.6202 0.2357 0.1980 -0.2393 0.0155 -0.0890
0.2950 0.4295 -0.4772 -0.3318 -0.0880 -0.3618
-0.4784 0.3112 -0.3427 -0.2697 0.2311 0.1248
0.2629 0.3668 -0.0962 -0.1104 -0.5778 0.4586
0.1937 0.4314 0.2660 -0.0783 0.4344 -0.0199
0.2279 0.3233 0.4762 0.0706 0.0093 -0.4843
-0.1725 0.2991 -0.3651 0.8274 -0.0094 -0.1772
-0.3316 0.0623 0.2617 0.0026 -0.6411 -0.2801
Check if $Im(B) \subset X_c$.
Background. Let $V$ and $W$ be two subspaces of $\mathbb{R}^n$, dim$(V) = k$, dim$(W) = m$. $V \subset W$ if there exists a matrix
$$ A = \begin{pmatrix} a_{11} & ... & a_{1k} \\ ... & ... & ... \\ a_{m1} & ... & a_{mk} \end{pmatrix} $$
such that $V = W A$. If there exists such matrix, it can be computed as follows: $A = \left(W^T W\right)^{-1} W^T V$. If the equility $V = W A$ holds for this matrix $A$, than $V \subset W$. So, we need that $W \left(W^T W\right)^{-1} W^T V - V$ be approximately the zero matrix.
For this purpose, we define the following function (in a separate file, called pcz_vecalg_subset.m)
function [ret] = pcz_vecalg_subset(V, W, tol)
ret = rank(W / (W'*W) * W' * V - V, tol) == 0;
Using this function we can check whether $Im(B) \subset X_c$.
pcz_info(pcz_vecalg_subset(orth(B),X_C,1e-10), 'Im(B) ⊂ X_C')
│ [ OK ] Im(B) ⊂ X_C
Check if $A X_c \subseteq X_c$.
pcz_info(pcz_vecalg_subset(A*X_C,X_C,1e-10), 'A X_C ⊆ X_C')
disp '$a$'
│ [ OK ] A X_C ⊆ X_C $a$