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
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
|
function []=evans(n,d,kmax)
//seuil maxi et mini (relatifs) de discretisation
// en espace
// Copyright INRIA
smax=0.01;smin=smax/3;
nptmax=500 //nbre maxi de pt de discretisation en k
//
//analyse de la liste d'appel
[lhs,rhs]=argn(0)
if rhs<=0 then //demo
s_mat=['xbasc();n=real(poly([0.1-%i 0.1+%i,-10],''s''));';
' d=real(poly([-1 -2 -%i %i],''s''));';
'evans(n,d,80);'];
write(%io(2),s_mat);execstr(s_mat);
return;
end
select type(n)
case 1 then
if rhs==2 then kmax=0,end
case 2 then
if rhs==2 then kmax=0,end
//-compat next case retained for list/tlist compatibility
case 15 then
if rhs==2 then kmax=d,else kmax=0,end
n1=n(1);
select n1(1)
case 'r' then [n,d]=n(2:3)
case 'lss' then n=ss2tf(n);[n,d]=n(2:3);n=clean(n);d=clean(d);
else error('transfer or state-space only')
end
case 16 then
if rhs==2 then kmax=d,else kmax=0,end
n1=n(1);
select n1(1)
case 'r' then [n,d]=n(2:3)
case 'lss' then n=ss2tf(n);[n,d]=n(2:3);n=clean(n);d=clean(d);
else error('transfer or state-space only')
end
else error('transfer or state-space only')
end
if prod(size(n))<>1 then
error('SISO system only'),
end
if kmax<=0 then
nm=mini([degree(n),degree(d)])
fact=norm(coeff(d),2)/norm(coeff(n),2)
kmax=round(500*fact),
end
//
//calcul de la discretisation en k et racines associees
nroots=roots(n);racines=roots(d);
if nroots==[] then
nrm=maxi([norm(racines,1),norm(roots(d+kmax*n),1)])
else
nrm=maxi([norm(racines,1),norm(nroots,1),norm(roots(d+kmax*n),1)])
end
md=degree(d)
//
ord=1:md;kk=0;nr=1;k=0;pas=0.99;fin='no';
while fin=='no' then
k=k+pas
r=roots(d+k*n);r=r(ord)
dist=maxi(abs(racines(:,nr)-r))/nrm
//
point='nok'
if dist <smax then //pas correct
point='ok'
else //pas trop grand ou ordre incorrect
// on cherche l'ordre qui minimise la distance
ix=1:md
ord1=[]
for ky=1:md
yy=r(ky)
mn=10*dist*nrm
for kx=1:md
if ix(kx)>0 then
if abs(yy-racines(kx,nr)) < mn then
mn=abs(yy-racines(kx,nr))
kmn=kx
end
end
end
ix(kmn)=0
ord1=[ord1 kmn]
end
r(ord1)=r
dist=maxi(abs(racines(:,nr)-r))/nrm
if dist <smax then
point='ok',
ord(ord1)=ord
else
k=k-pas,pas=pas/2.5
end
end
if dist<smin then
//pas trop petit
pas=2*pas;
end
if point=='ok' then
racines=[racines,r];kk=[kk,k];nr=nr+1
if k>kmax then fin='kmax',end
if nr>nptmax then fin='nptmax',end
end
end
//taille du cadre graphique
x1 =[nroots;matrix(racines,md*nr,1)];
xmin=mini(real(x1));xmax=maxi(real(x1))
ymin=mini(imag(x1));ymax=maxi(imag(x1))
dx=abs(xmax-xmin)*0.05
dy=abs(ymax-ymin)*0.05
rect=[xmin-dx;ymin-dy;xmax+dx;ymax+dy];
dx=maxi(abs(xmax-xmin),abs(ymax-ymin));
//trace des lieux des zeros
xx=xget("mark")
xset("mark",xx(1),xx(1)+3);
if nroots<>[] then
plot2d(real(nroots),imag(nroots),[-5,4],"151",...
'open loop zeroes',rect);
strf="100"
else
strf="151"
end
//trace des lieu des poles en boucle ouverte
if racines<>[] then
plot2d(real(racines(:,1)),imag(racines(:,1)),[-2,5],strf,...
'open loop poles',rect);
strf='100';else strf='151';
end
// trace des branches asymptotiques
plot2d(0,0,[1,2],strf,'asymptotic directions');
xtitle('Evans root locus','Real axis','Imag. axis');
xset("clipgrf");
m=degree(n);q=md-m
if q>0 then
la=0:q-1;
so=(sum(racines(:,1))-sum(nroots))/q
i1=real(so);i2=imag(so);
if prod(size(la))<>1 then
ang1=%pi/q*(ones(la)+2*la)
x1=dx*cos(ang1),y1=dx*sin(ang1)
// ang2=%pi/q*2*la
// x2=dx*cos(ang2),y2=dx*sin(ang2)
else
x1=0,y1=0,//x2=0,y2=0,
end
if md==2,
if coeff(d,md)<0 then
x1=0*ones(2),y1=0*ones(2)
// x2=0*ones(2),y2=0*ones(2),
end,
end;
if maxi(k)>0 then
for i=1:q,xsegs([i1,x1(i)+i1],[i2,y1(i)+i2]),end,
end
// if mini(k)<0 then
// for i=1:q,xsegs([i1,x2(i)+i1],[i2,y2(i)+i2]),end,
// end
end;
xclip();
//lieu de evans
[n1,n2]=size(racines);
plot2d(real(racines)',imag(racines)',3*ones(1,n2),"100",...
'closed-loop poles loci');
if fin=='nptmax' then
write(%io(2),'evans : too many points required')
end
xset("mark",xx(1),xx(2));
// gain corresponding to a selected point of the locus
//[l1,l2]=min(abs(racines(:)-selected));
//col=int(l2/n1)+1;//row=modulo(l2,n1); racines(row,col) <-> seleceted
//gain=kk(col);
|
