这是程序代码:
E=10;
L1=15e-6;
C1=100e-6;
L2=15e-6;
C2=100e-6;
R=20;
r=10e-3;
Iref1=1.3;
Ton=0.3e-6;
Iref2=1.5;
T2=0.88e-6;
Pcpl=8;
I=[1,0,0,0;0,1,0,0;0,0,1,0;0,0,0,1];
syms Vo1 Vo2 iL1 iL2
for i=1:100
Pcpl=8+(i-1)*0.02;%变化的参数
Id=iL1+Ton*(E-Vo1)/L1;
iL2=Pcpl/Vo1;
d1=Ton/(Ton+(Id-Iref1)*L1/Vo1);
d2=(Iref2-iL2)*L2/((Vo1-Vo2)*T2);
%d1(d1>=1)=0.99;d1(d1<=0)=0.01;
%d2(d2>=1)=0.99;d2(d2<=0)=0.01;
% A1=[0,-1/L1,0,0;r+1/C1,0,0,0;0,0,0,-1/L2;0,0,1/C2,-1/(R*C2)];
A1=[0,-1/L1,0,0;r+1/C1,0,(-r-1/C1),0;0,1/L2,0,-1/L2;0,0,1/C2,-1/(R*C2)];←←←←
B1=[d1/L1,0,0,0]';
A2=[0,-1/L1,0,0;r+1/C1,0,0,0;0,0,0,-1/L2;0,0,1/C2,-1/(R*C2)];
B2=[d1/L1,0,0,0]';
Xn=[iL1;Vo1;iL2;Vo2];
X=expm(A2*(1-d2)*T2)*expm(A1*d2*T2)*Xn+expm(A2*(1-d2)*T2)*(expm(A1*d2*T2)-I)*inv(A1)*B1*E+(expm(A2*(1-d2)*T2)-I)*inv(A2)*B2*E;% 非线性函数系统
f=Xn-X;
fixed_point=fsolve(f,[1 1 1 1]);%平衡点 ←←←←
jacobian_mat=jacobian(X,[iL1;Vo1;iL2;Vo2]);%雅可比矩阵
%n=length(fixed_point.iL1);%n为平衡点组数,即n组平衡点[x1 x2]
n=size(fixed_point,1);
for j=1:n
fixed(j).jacobian=subs(jacobian_mat,Xn,{fixed_point.Xn(j)});%每个平衡点的雅可比矩阵
fixed(j).eig=eig(fixed(j).jacobian);%平衡点的特征值
eig_max=max(double(real(fixed(j).eig)));
if eig_max<=0
plot(Pcpl,double(fixed_point.Xn(j)),'.')
hold on
else
plot(Pcpl,double(fixed_point.Xn(j)),'+')%鞍点用+标出
hold on
end;
end;
end
经过试验发现是矩阵A1无法带指数e计算,一直正忙,后来我尝试换成上面带%位置的A1矩阵又可以了,请问应该如何解决?
另外这是一个求非线性方程平衡点数值解、雅可比矩阵特征值的程序,我用solve就会警告算不出结果,用如上的fsolve就报错:如果 FUN 为 MATLAB 对象,则它必须具有 feval 方法。请问这是问题出在哪?
求各位帮忙分析分析,感谢!