hello everyone I have bellow code to solve with ode solvers but souloution goes to infinity.I have checked the model and I'm sure about its correctness.thanks for every help

function out=vbrcont2(t,s)

theta=s(1);
omegha=s(2);
lambdafd=s(3);
lambdakd1=s(4);
lambdakq1=s(5);
lambdakq2=s(6);
ia=s(7);
ib=s(8);
ic=s(9);
alpha=theta;%phase a angle
beta=theta-2*pi/3;% phase b angle
ghama=theta+2*pi/3;%phase c angle

VL=26e3;
p=2;
we=2*pi*60;
rs=0.00234;                                                                 
Xls=0.1538;                                                                 
Lls=Xls/we;
Xq=1.457;
Xd=1.457;
Xmd=Xd-Xls;
Lmd=Xmd/we;
Xmq=Xq-Xls;
Lmq=Xmq/we;
rkq1=0.00144;
Xplkq1=0.6578;
Llkq1=Xplkq1/we;
rkq2=0.00681;
Xplkq2=0.07602;
Llkq2=Xplkq2/we;
rfd=0.00075;
Xplfd=0.1145;
Llfd=Xplfd/we;
rkd1=0.01080;
Xplkd=0.06577;
Llkd1=Xplkd/we;
J=0.0658e6;

vfd=sqrt(2/3)*VL*(rfd/Xmd);
Tm=0;

%/////////////////////////////////////////////////////////
park_T=(2/3)*[cos(alpha) cos(beta) cos(ghama);sin(alpha) sin(beta) sin(ghama);
    0.5 0.5 0.5 ];
iss=park_T*[ia;ib;ic];
iqs=iss(1,1);
ids=iss(2,1);
Lmqdp=(((1/Lmq)+(1/Llkq1)+(1/Llkq2))^-1);%dp means double prime(subtransiant) 
Lmddp=(((1/Lmd)+(1/Llfd)+(1/Llkd1))^-1);%dp means double prime(subtransiant)
La=(Lmqdp+Lmddp)/3;
Lb=(Lmddp-Lmqdp)/3;
lambdamq=Lmqdp*((lambdakq1/Llkq1)+(lambdakq2/Llkq2)+iqs);
lambdamd=Lmddp*((lambdafd/Llfd)+(lambdakd1/Llkd1)+ids);
Te=(3*p/4)*(lambdamd*iqs-lambdamq*ids);
lambdaqdp=Lmqdp*((lambdakq1/Llkq1)+(lambdakq2/Llkq2));
lambdaddp=Lmddp*((lambdafd/Llfd)+(lambdakd1/Llkd1));
vqdp=(omegha*lambdaddp)+((Lmqdp*rkq1/(Llkq1^2))*(lambdaqdp-lambdakq1))+((Lmqdp*rkq2/(Llkq2^2))*(lambdaqdp-lambdakq2))+(((rkq1/(Llkq1^2))+(rkq2/(Llkq2^2)))*((Lmqdp^2)*iqs));
vddp=(-omegha*lambdaqdp)+((Lmddp*rkd1/(Llkd1^2))*(lambdaddp-lambdakd1))+((Lmddp/Llfd)*vfd)+((Lmddp*rfd/(Llfd^2))*(lambdaddp-lambdafd))+(((rfd/(Llfd^2))+(rkd1/(Llkd1^2)))*((Lmddp^2)*ids));
Inv_park_T=[cos(alpha) sin(alpha) 1;cos(beta) sin(beta) 1;cos(ghama) sin(ghama) 1];
vdp=Inv_park_T*[vqdp;vddp;0];
vadp=vdp(1,1);
vbdp=vdp(2,1);
vcdp=vdp(3,1);
va=VL*sqrt(2/3)*cos(we*t);
vb=VL*sqrt(2/3)*cos(we*t-(2*pi/3));
vc=VL*sqrt(2/3)*cos(we*t+(2*pi/3));
Labcdp=1*[Lls+La-Lb*cos(2*theta),-(La/2)-Lb*cos(2*theta-(2*pi/3)),-(La/2)-Lb*cos(2*theta+(2*pi/3));
    -(La/2)-Lb*cos(2*theta-(2*pi/3)),Lls+La-Lb*cos(2*theta-(4*pi/3)),-(La/2)-Lb*cos(2*theta);
    -(La/2)-Lb*cos(2*theta+(2*pi/3)),-(La/2)-Lb*cos(2*theta),Lls+La-Lb*cos(2*theta+(4*pi/3))];
% inv_Labcdp=inv(Labcdp);

% Labcdp*inv_Labcdp
Y=(1/1)*[-va-rs*ia-2*Lb*omegha*sin(2*theta)*ia-2*Lb*omegha*sin(2*theta-(2*pi/3))*ib-2*Lb*omegha*sin(2*theta+(2*pi/3))*ic+vadp
    ;-vb-rs*ib-2*Lb*omegha*sin(2*theta-(2*pi/3))*ia-2*Lb*omegha*sin(2*theta-(4*pi/3))*ib-2*Lb*omegha*sin(2*theta)*ic+vbdp
    ;-vc-rs*ic-2*Lb*omegha*sin(2*theta+(2*pi/3))*ia-2*Lb*omegha*sin(2*theta)*ib-2*Lb*omegha*sin(2*theta+(4*pi/3))*ic+vcdp];
% det(Labcdp);
di=Labcdp\Y;

iadot=di(1,1);
ibdot=di(2,1);
icdot=di(3,1);
out=[omegha;-(p/(2*J))*(Te-Tm);-(rfd/Llfd)*(lambdafd-lambdamd)+vfd;-(rkd1/Llkd1)*(lambdakd1-lambdamd)
    ;-(rkq1/Llkq1)*(lambdakq1-lambdamq);-(rkq2/Llkq2)*(lambdakq2-lambdamq);iadot;ibdot;icdot];
end

[Walter Roberson]

We have no hope of examining the correctness of this code without seeing the equations that are to be solved.

As a practical matter, I am also not going to attempt to debug code that long that does not use whitespace to make the expressions clearer. Spaces do not cost anything to put in, but they make it much easier for other people to read. I recommend a space before and after each operator, so for example

((p+q)./r-s*t)

would become

((p + q) ./ r - s * t)

[behnam ojaghloo]

Hi dear walter.This is what I'm going to solve(pdf named vbr). and here is parameters needed.

rs=0.00234;                                                                
Xls=0.1538;                                                                
Lls=Xls/we;
Xq=1.457;
Xd=1.457;
Xmd=Xd-Xls;
Lmd=Xmd/we;
Xmq=Xq-Xls;
Lmq=Xmq/we;
rkq1=0.00144;
Xplkq1=0.6578;
Llkq1=Xplkq1/we;
rkq2=0.00681;
Xplkq2=0.07602;
Llkq2=Xplkq2/we;
rfd=0.00075;
Xplfd=0.1145;
Llfd=Xplfd/we;
rkd1=0.01080;
Xplkd=0.06577;
Llkd1=Xplkd/we;

for equations 29 & 30 is 2

Answer 1

[Jan]

I'm convinced, that it is impossible to debug this code. Only a genious can keep the overview in this pile of expressions. Therefore I recommend to split the formula into parts, which results can be validated. There is simply no other chance to debug the code.

Expression which appeare repeatedly should be stored in a constant like theta-(2*pi/3) . Other constants like VL=26e3 are defined, but a 26000 appears in the code also. Expressions like "(1/1)*" and "1*" should be removed. Typing matrices in their real 2D format increases the readability also. I like the usage of parenthesis, even if they are not required, but in these line the outer ones are only optical junk:

La=((Lmqdp+Lmddp)/3);
va=(26000*sqrt(2/3))*cos(we*t);

The multiplication with the inverse of a matrix is not a smart calculation from a numerical point of view. Try to calculated the transposed result using the \ operator.

The debugging will be much easier, if the code is substantially simplified. And then a comparison with the mathematical formula will be easier - or better: possible for human.