diff options
Diffstat (limited to 'modules/cacsd/macros/copfac.sci')
| -rwxr-xr-x | modules/cacsd/macros/copfac.sci | 101 |
1 files changed, 101 insertions, 0 deletions
diff --git a/modules/cacsd/macros/copfac.sci b/modules/cacsd/macros/copfac.sci new file mode 100755 index 000000000..23536ed7d --- /dev/null +++ b/modules/cacsd/macros/copfac.sci @@ -0,0 +1,101 @@ +// 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.1-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), + if rhs<1 then + error(msprintf(gettext("%s: Wrong number of input arguments: At least %d expected.\n"),"copfac",1)) + end + select typeof(g) + case "rational" then + g=tf2ss(g) + case "state-space" then + else + error(msprintf(gettext("%s: Wrong type for input argument #%d: Linear state space or a transfer function expected.\n"),"copfac",1)) + end + if g.dt<>"c" then + error(msprintf(gettext("%s: Wrong type for argument #%d: In continuous time expected.\n"),"copfac",1)) + end + + [r,p,t]=size(g); + [a,b,c,d]=abcd(g), + [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 + itol=2 + tol=polf,polc=-ones(1,n1),polf=-ones(1,n2), + case 3 then + tol=1000*%eps, + else + itol=4 + end, + if type(tol)<>1|size(tol,"*")<>1|~isreal(tol)|tol<=0 then + error(msprintf(gettext("%s: Wrong value for input argument #%d: Must be a positive scalar.\n"),"copfac",itol)) + end + if type(polf)<>1 then + error(msprintf(gettext("%s: Wrong type for input argument #%d: Vector of floats expected.\n"),"copfac",2)) + end + if size(polf,"*")<>n2 then + error(msprintf(gettext( "%s: Wrong size for argument #%d: %d expected.\n"),"copfac",2,n2)) + end + if type(polc)<>1 then + error(msprintf(gettext("%s: Wrong type for input argument #%d: Vector of floats expected.\n"),"copfac",3)) + end + if size(polc,"*")<>n1 then + error(msprintf(gettext( "%s: Wrong size for argument #%d: %d expected.\n"),"copfac",3,n1)) + end + + + //-------------------- + if n1<>t then + w1=u1(:,n1+1:t),a1=w1'*a*w1, + no=norm(a1), + if max(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 max(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 |