0

我们正在研究有一个求解器的项目,该求解器使用 Silvester 矩阵来求解丢番图方程,但结果不是我们预期的,似乎 Silvester 矩阵代码存在问题。

主要代码:

% The system transfer function
C =[0 0 1]; % numerator ( descending powers )
D =[1 0.2 1]; % denominator ( descending powers )
% ... The close loop characteristics :
% Settling time ts ~ 4s
ts =4;
% tau = ts /(2.5~3)
tau = ts /2.5;
% Able to completely suppress step disturbances .
% => controller order (n) = 2
n =2;
% Characteristic Polynomial Order (q) = 4
q =4;
 
gama =[2.5 2 2];
% deriving the characteristic polynomial
P= charactpoly (tau , gama ,q);
% Computing the Sylvester Matrix
S= sylvester4cdm (C ,D ,n ,1) ;

X=S\P. ’;

% The numerator and denominator
A = [(X(1:2:end-1)).' , 0]
B = [(X(2:2:end)).' , X(end)]

% The pre -filter (just a scaling constant)
F = polyval(P,0)/polyval(C,0)

charactpoly 功能代码:

function P= charactpoly (tau , gama , np )
 %
 % np - the characteristic polynomial order
 % tau - the equivalent time constant
 % gama - the vector of stability values
 % P - the characteristic polynomial vectorin decreasing
 % order . i.e. P (1) *s^ np + ... + P(end -1)*s + P( end )
 %
 % e.g. charactpoly (1 ,[2.5 2 2 2.5] ,5)

 P= zeros (1 , np +1) ;
 P(np +1) =1; P(np)= tau ;
 for i =2: np
 j =1:i -1;
 % The characteristic coefficients vector
 P(np -i +1) = tau ^i* prod (1./( gama (j)) .^(i -j));
 end

Sylvester4cdm 功能:

function S = cdmsylvester(C,D,n,k)
 
 % n : Controller order
 % k : Number of known A(m) coefficients
 % C : Transfer function numerator coefficients
 % D : Transfer function denominator coefficients
 % S : The Sylvester matrix (m+n+1 x 2n-k+2)

 m = length(D)-1; % System order
 nl = m+n-k+1;
 nc = (2*n)-(2*k)+2;
 
 Df =[D zeros(1,n-k)];
 
 Cf =[C zeros(1,n-k)];
 
 A = toeplitz(Df ,[Df(1)  zeros(1,n-k)]);
 
 B = toeplitz(Cf ,[Cf(1)  zeros(1,n-k)]);
 
 S1 = reshape([A;B],nl,nc);
 
 % Adding the necessary zero rows to matrix
 
 S2 = [zeros(k,nc);S1];
 
 % Adding the columns due only to B(m)
 
 Cf = fliplr(C);
 
 S3 = toeplitz ([Cf  zeros(1,n)],[Cf(1)  zeros(1,k-1)]);
 
 if k>0,S=[S3 S2]; else S=S1; end;

运行主代码后的最终答案应该是:

A(s) = 0.0524s^2 + 0.3172s & B(s) = 0.908s^2 + 1.283s + 1 & F(s) = 1

如果您能帮助我们解决这个问题,我们将不胜感激。

4

0 回答 0