1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
|
function [X,dims,J,Y,k,Z]=cainv(Sl,alfa,beta,flag)
//CA invariant subspace: dual of abinv
//Finds orthogonal bases X and Y and output injection J
//such that the abcd matrices of Sl in bases (X,Y) are displayed as:
//
// [A11,*,*,*,*,*] [*]
// [0,A22,*,*,*,*] [*]
// X'*(A+J*C)*X = [0,0,A33,*,*,*] X'*(B+J*D) = [*]
// [0,0,0,A44,*,*] [0]
// [0,0,0,0,A55,*] [0]
// [0,0,0,0,0,A66] [0]
//
// Y*C*X = [0,0,C13,*,*,*] Y*D = [*]
// [0,0,0,0,0,C26] [0]
//
// dims=[nd1,nu1,dimS,dimSg,dimN] defines the partition of A matrix and k the partition
// of [C,D] matrix.
// eigenvalues of A11 (nd1 x nd1) are unstable
// eigenvalues of A22 (nu1-nd1 x nu1-nd1) are stable
// pair (A33, C13) (dimS-nu1 x dimS-nu1, k x dimS-nu1) is observable, eigenvalues of A33 set to alfa
// A44 (dimSg-dimS x dimSg-dimS) is unstable
// A55 (dimN-dimSg,dimN-dimSg) is stable
// pair (A66,C26) (nx-dimN x nx-dimN, ) is observable, eigenvalues of A66 set to beta.
//
// dimS first columns of X span S= smallest (C,A) invariant
// subspace which contains Im(B).
// dimSg first columns of X span Sg
// dimN first columns of X span N=S+V
// dimS=0 iff B(ker(D))=0
//
// DDEP: dot(x)=A x + Bu + Gd
// y= Cx (observation)
// z= Hx (z=variable to be estimated, d=disturbance)
// Find: dot(w) = Fw + Ey + Ru such that
// zhat = Mw + Ny
// z-Hx goes to zero at infinity
// Solution exists iff Ker H contains Sg(A,C,G) inter KerC
//i.e. H is such that:
// For any W which makes a column compression of [Xp(1:dimSg,:);C]
// with Xp=X' and [X,dims,J,Y,k,Z]=cainv(syslin('c',A,G,C));
// [Xp(1:dimSg,:);C]*W = [0 | *] one has
// H*W = [0 | *] (with at least as many columns as above).
// Copyright INRIA
[LHS,RHS]=argn(0);
if RHS==1 then alfa=-1;beta=-1;flag='ge';end
if RHS==2 then beta=alfa;flag='ge';end
if RHS==3 then flag='ge';end
if RHS==4 then
if type(flag)~=10 then error('abinv: flag must be a string');end
end
[X,ddims,F,U,k,Z]=abinv(Sl',beta,alfa,flag);
[nx,nx]=size(X);
select flag
case 'ge'
nr=ddims(1);nvg=ddims(2);nv=ddims(3);noc=ddims(4);nos=ddims(5);
nd1=nx-nos;nu1=nx-noc;dimS=nx-nv;dimSg=nx-nvg;dimN=nx-nr;
n6=1+ddims(5):nx;
n5=1+ddims(4):ddims(5);
n4=1+ddims(3):ddims(4);
n3=1+ddims(2):ddims(3);
n2=1+ddims(1):ddims(2);
n1=1:ddims(1);
//nr=1:nr;nzs=nr+1:nr+nvg;nzi=nvg+1:nv;
X=[X(:,n6),X(:,n5),X(:,n4),X(:,n3),X(:,n2),X(:,n1)];
J=F';Z=Z';Y=U';Y=[Y(k+1:$,:);Y(1:k,:)];
dims=[nd1,nu1,dimS,dimSg,dimN];
return;
case 'st'
dims=nx-ddims;dims=dims($:-1:1);
n5=1+ddims(4):nx;
n4=1+ddims(3):ddims(4);
n3=1+ddims(2):ddims(3);
n2=1+ddims(1):ddims(2);
n1=1:ddims(1);
X=[X(:,n5),X(:,n4),X(:,n3),X(:,n2),X(:,n1)];
J=F';Z=Z';Y=U';Y=[Y(k+1:$,:);Y(1:k,:)];
return;
case 'pp'
dims=nx-ddims;dims=dims($:-1:1);
n4=1+ddims(3):nx;
n3=1+ddims(2):ddims(3);
n2=1+ddims(1):ddims(2);
n1=1:ddims(1);
X=[X(:,n4),X(:,n3),X(:,n2),X(:,n1)];
J=F';Z=Z';Y=U';Y=[Y(k+1:$,:);Y(1:k,:)];
return;
end
|
