1// Scilab ( http://www.scilab.org/ ) - This file is part of Scilab 2// Copyright (C) 1996 - 2016 - INRIA - Serge Steer 3// Copyright (C) 2012 - 2016 - Scilab Enterprises 4// 5// This file is hereby licensed under the terms of the GNU GPL v2.0, 6// pursuant to article 5.3.4 of the CeCILL v.2.1. 7// This file was originally licensed under the terms of the CeCILL v2.1, 8// and continues to be available under such terms. 9// For more information, see the COPYING file which you should have received 10// along with this program. 11 12function [kp,s]=krac2(sys) 13 //The denominator of the closed loop system is den(s)+K*num(s). So the 14 // the closed loops poles verify K(s)=-den(s)/num(s) 15 //The real axis breakaway points occurs at the extrema of the K(s) 16 // so at the point where K'=dK/ds = 0 17 // K'=-(den'*num-den*num')/num^2 18 // K'= 0 --> den'*num-den*num'=0 19 // http://www.scribd.com/doc/21374148/An-Introduction-to-Control-Systems 20 if and(typeof(sys)<>["state-space" "rational" "zpk"]) then 21 ierr=execstr("[kp,s]=%"+typeof(sys,"overload")+"_krac2(sys)","errcatch") 22 if ierr<>0 then 23 error(msprintf(_("%s: Wrong type for input argument #%d: Linear dynamical system expected.\n"),"krac2",1)) 24 end 25 return 26 end 27 28 if size(sys,"*")<>1 then 29 error(msprintf(gettext("%s: Wrong size for input argument #%d: Single input, single output system expected.\n"),"krac2",1)) 30 end 31 32 if typeof(sys)=="state-space" then 33 sys=ss2tf(sys) 34 elseif typeof(sys)=="zpk" then 35 sys=zpk2tf(sys) 36 end 37 num=sys.num; 38 den=sys.den; 39 40 s=roots(derivat(num)*den-derivat(den)*num,"e") 41 //collect the real roots only 42 i=find(abs(imag(s))<=10*%eps) 43 if i==[] then kp=[],s=[];return,end 44 s=real(s(i))' 45 s=s(horner(num,s)<>0); 46 kp=-real(freq(den,num,real(s))); 47 i=find(kp>=0); 48 kp=kp(i) 49 s=s(i) 50endfunction 51 52