1// Scilab ( http://www.scilab.org/ ) - This file is part of Scilab 2// Copyright (C) INRIA - 3// 4// Copyright (C) 2012 - 2016 - Scilab Enterprises 5// 6// This file is hereby licensed under the terms of the GNU GPL v2.0, 7// pursuant to article 5.3.4 of the CeCILL v.2.1. 8// This file was originally licensed under the terms of the CeCILL v2.1, 9// and continues to be available under such terms. 10// For more information, see the COPYING file which you should have received 11// along with this program. 12 13function [n]=linf(g,eps,tol) 14 //linf(g,[eps],[tol]) L_infinity norm 15 // n=sup [sigmax(g(jw)] (sigmax largest singular value). 16 // w 17 //-- g is a syslin system. 18 //-- eps is error tolerance on n. 19 //-- tol threshold for imaginary axis poles. 20 // See also: h_norm 21 //! 22 if type(g)==1,if norm(g)==0,n=0,return,end,end, 23 24 if and(typeof(g)<>["rational","state-space"]) then 25 error(msprintf(gettext("%s: Wrong type for input argument #%d: Linear state space or a transfer function expected.\n"),"linf",1)) 26 end 27 if g.dt<>"c"&g.dt<>[] then 28 error(msprintf(gettext("%s: Wrong value for input argument #%d: Continuous time system expected.\n"),"linf",1)) 29 end 30 g.dt="c" 31 if typeof(g)=="rational" then g=tf2ss(g),end 32 33 [lhs,rhs]=argn(0), 34 select rhs, 35 case 1 then eps=1e-7,tol=1000*%eps, 36 case 2 then tol=1000*%eps, 37 end, 38 [a,b,c,d]=g(2:5),[t,t]=size(a), 39 p=ctr_gram(g),q=obs_gram(g); 40 //Algorithm: 41 //---------- 42 //1. min , max. 43 //---------------------------------- 44 [slp,slm]=dtsi(g), 45 if slp==0 then pp=0,qq=0,tp=1, 46 pm=p,qm=q,tm=t, 47 else 48 if slm==0 then pm=0,qm=0,tm=1, 49 pp=p,qq=q,tp=t, 50 else 51 [tp,tp]=size(slp(2)),[tm,tm]=size(slm(2)), 52 pp=ctr_gram(slp),qq=obs_gram(slp), 53 pm=ctr_gram(slm),qm=obs_gram(slm), 54 end, 55 end, 56 hsvp=sqrt(spec(pp*qq)),hsvp=gsort(real(hsvp)), 57 hsvm=sqrt(spec(pm*qm)),hsvm=gsort(real(hsvm)), 58 gl=max([norm(d),hsvp(tp),hsvm(tm)]), 59 gu=norm(d)+2*(sum(hsvp)+sum(hsvm)), 60 //2. binary search 61 //---------------------- 62 while gu-gl>2*eps*gl, 63 x=(gl+gu)/2, 64 r=d'*d-(x*x)*eye(),s=d*d'-(x*x)*eye(), 65 mx=[a-b/r*d'*c, -x*b/r*b'; .. 66 x*c'/s*c, -a'+c'*d/r*b'], 67 mp=abs(real(spec(mx))),mp=min(mp), 68 if mp>tol then gu=x, else gl=x, end, 69 end; 70 n=(gu+gl)/2 71endfunction 72