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
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
|
function [h,num,den]=ss2tf(sl,rmax)
// State-space to transfer function.
// Syntax:
// h=ss2tf(sl)
// h=ss2tf(sl,'b')
// h=ss2tf(sl,rmax)
//
// sl : linear system (syslin list)
// h : transfer matrix
// rmax : optional parameter controlling the conditioning in
// block diagonalization method is used (If 'b' is entered
// a default value is used)
// Method: By default, one uses characteristic polynomial and
// det(A+Eij)=det(A)+C(i,j) where C is the adjugate matrix of A
// Other method used : block-diagonalization (generally
// this gives a less accurate result).
//
//!
// Copyright INRIA
if type(sl)==1|type(sl)==2 then
h=sl
return
end
flag=sl(1);
if (type(sl)<>16)|flag(1)<>'lss' then
error('First argument must be in state-space form')
end
if sl(3)==[] then h=sl(5);num=sl(5);den=eye(sl(5));return;end
if sl(4)==[] then h=sl(5);num=sl(5);den=eye(sl(5));return;end
if size(sl(2),'*')==0 then
h=sl(5)
return
end
//1
domaine=sl(7)
if type(domaine)==1 then var='z';end
if domaine=='c' then var='s';end;
if domaine=='d' then var='z';end;
if domaine==[] then
var='s';
if type(sl(5))==2 then var=varn(sl(5));end
end
//
[lhs,rhs]=argn(0);
meth='p'
if rhs==2 then
if type(rmax)==10 then
meth=part(rmax,1),
else
meth='b'
end
end
//
select meth
case 'b' //
a=sl(2)
[n1,n1]=size(a)
z=poly(0,var);
if rhs==1 then
[a,x,bs]=bdiag(a)
else
[a,x,bs]=bdiag(a,rmax)
end
k=1;m=[];v=ones(1,n1);den=1;
for n=bs';k1=k:k-1+n;
// leverrier
h=z*eye(n,n)-a(k1,k1)
f=eye(n,n)
for kl=1:n-1,
b=h*f,
d=-sum(diag(b))/kl,
f=b+eye()*d,
end;
d=sum(diag(h*f))/n
//
den=den*d;
l=[1:k-1,k+n:n1] ,
if l<>[] then v(l)=v(l)*d,end
m=[m,x(:,k1)*f];
k=k+n;
end;
if lhs==3 then h=sl(5),num=sl(4)*m*diag(v)*(x\sl(3));return;end
m=sl(4)*m*diag(v)*(x \ sl(3))+sl(5)*den;
[m1,n1]=size(m);[m,den]=simp(m,den*ones(m1,n1))
h=syslin(domaine,m,den)
case 'p' then
Den=real(poly(sl(2),var))
na=degree(Den);den=[];
[m,n]=size(sl(5))
c=sl(4)
for l=1:m
[m,i]=maxi(abs(c(l,:)));
if m<>0 then
ci=c(l,i)
t=eye(na,na)*ci;t(i,:)=[-c(l,1:i-1), 1, -c(l,i+1:na)]
al=sl(2)*t;
t=eye(na,na)/ci;t(i,:)=[c(l,1:i-1)/ci, 1, c(l,i+1:na)/ci]
al=t*al;ai=al(:,i),
b=t*sl(3)
for k=1:n
al(:,i)=ai+b(:,k);
[nlk,dlk]=simp(real(poly(al,var)),Den)
den(l,k)=dlk;
num(l,k)=-(nlk-dlk)*ci
end
else
num(l,1:n)=0*ones(1,n);
den(l,1:n)=ones(1,n);
end
end
if lhs==3 then h=sl(5);return;end
w=num./den+sl(5);
if type(w)==1 then h=w;return;end //degenerate case
h=syslin(domaine,w);
end
|
