我们正在研究有一个求解器的项目,该求解器使用 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
如果您能帮助我们解决这个问题,我们将不胜感激。