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
|
// Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
// Copyright (C) INRIA -
//
// This file must be used under the terms of the CeCILL.
// This source file is licensed as described in the file COPYING, which
// you should have received as part of this distribution. The terms
// are also available at
// http://www.cecill.info/licences/Licence_CeCILL_V2-en.txt
function [n,m,xt,yt]=copfac(g,polf,polc,tol)
//[n,m,xt,yt]=copfac(G,[polf,polc,[tol]]) returns a right coprime
//factorization of g :
// g = n*m^-1 where n and m are stable, proper and right coprime.
// (ie. [n m] left-invertible with stability)
// xt and yt satisfy:
// [xt -yt].[m n]' = eye (Bezout identity)
// G is assumed stabilizable and detectable.
//-- G is is linear system (syslin list)
//-- polf are polc respectively the poles of xt and yt
// and the poles of n and m .
// These are optional arguments with defautl values -1.
//-- tol is a threshold for detecting stable poles.
//!
[lhs,rhs]=argn(0),
select typeof(g)
case "rational" then
g=tf2ss(g)
case "state-space" then
else
error(msprintf(_("%s: Wrong type for input argument #%d: %s data structure expected.\n"),"copfac",1,"syslin"));
end
[r,p,t]=size(g);
[a,b,c,d]=g(2:5),
[n1,u1]=contr(a,b),[n2,u2]=contr(a',c'),
select rhs,
case 1 then
polc=-ones(1,n1),polf=-ones(1,n2),tol=1000*%eps,
case 2 then
tol=polf,polc=-ones(1,n1),polf=-ones(1,n2),
case 3 then tol=1000*%eps,
end,
//--------------------
if n1<>t then
w1=u1(:,n1+1:t),a1=w1'*a*w1,
no=norm(a1),
if maxi(real(spec(a1)))>no*tol then
error(msprintf(gettext("%s: Wrong value for input argument #%d: Stabilizable system expected.\n"),"copfac",1)),
end,
end,
//------------------
if n2<>t then
w2=u2(:,n2+1:t),a2=w2'*a*w2,
no=norm(a2),
if maxi(real(spec(a2)))>no*tol then
error(msprintf(gettext("%s: Wrong value for input argument #%d: Detectable system expected.\n"),"copfac",1)),
end,
end,
//-----------------------------
v1=u1(:,1:n1),a1=v1'*a*v1,b1=v1'*b,
k1=ppol(a1,b1,polc)*v1',
ak1=a-b*k1,ck1=c-d*k1,
v2=u2(:,1:n2),a2=v2'*a*v2,c2=c*v2,
k2=ppol(a2',c2',polf),
k2=v2*k2',
ak2=a-k2*c,bk2=b-k2*d,
//-------------------
m=syslin('c',ak1,b,-k1,eye(p,p)),
n=syslin('c',ak1,b,ck1,d),
xt=syslin('c',ak2,-bk2,-k1,eye(p,p)),
yt=syslin('c',ak2,k2,-k1),
endfunction
|
