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
|
// Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
// Copyright (C) INRIA - Serge Steer
//
// 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 [K,R]=kpure(sl,eps)
//if sl is a transfert function N(S)/D(s) kpure looks for K producing
//pure imaginary roots for
// P(s)=D(s)+K*N(s)
//There is a pair of pure imaginary poles if and only if
// P(%i*q)=0 (1)
// and
// P(-%i*q)=0 (2)
// because N and D are polynomials with real coefficients.
//Author: Serge Steer, INRIA
y=[];R=[];
msg=_("%s: Wrong type for input argument #%d: Linear state space or a transfer function expected.\n")
if argn(2)==1 then eps=1e-6,end
if size(eps,"*")==2 then eps=eps(2),end //compatibility
select typeof(sl)
case "rational" then
if size(sl.num,"*") <> 1 then
error(msprintf(msg,"kpure",1))
end
case "state-space" then
if size(sl.D,"*") <> 1 then
error(msprintf(msg,"kpure",1))
end
sl=ss2tf(sl)
else
error(msprintf(msg,"kpure",1))
end
//(1) give K(s)=-D(s)/N(s)
s=poly(0,varn(sl.den))
K=-sl.den/sl.num;
// replace K by the previous value in (2) and find the roots
s=roots(numer(horner(sl.den,-s)+K*horner(sl.num,-s)),"e");
//retain pure imaginary roots
s=imag(s(abs(real(s))<eps));
R=(s(s>0).'*%i);
//find the K(s) values K(s)=-D(s)/N(s)
K=-real(freq(sl.den,sl.num,R))
endfunction
|
