Coupled differential equation system

Hi

I tried to numerically solve a coupled differential equation system like this:

v=r*i1+L11*di1/dt+L12*di2/dt+L13*di3/dt;
v=r*i2+L21*di1/dt+L22*di2/dt+L23*di3/dt;
v=r*i3+L31*di1/dt+L32*di2/dt+L33*di3/dt;

[Image Analyst]

And did you succeed? If not read this link?

[siroos jalilyan]

No how i can solve with ode45 or other solver?

Answer 1

[Torsten]

Write your system as

L11*di1/dt+L12*di2/dt+L13*di3/dt = v-r*i1
L21*di1/dt+L22*di2/dt+L23*di3/dt = v-r*i2
L31*di1/dt+L32*di2/dt+L33*di3/dt = v-r*i3

and use the mass matrix option of the ode-solvers.

Best wishes

Torsten.

[siroos jalilyan]

thank you but i dont know how to solve with mass matrix option do you have any example that i use?

[Torsten]

L11=...;
L12=...;
L13=...;
L21=...;
L22=...;
L23=...;
L31=...;
L32=...;
L33=...;
L=[L11 L12 L13;L21 L22 L23;L31 L32 L33];
v=...;
r=...;
y0=...;
tspan=...;
fun=@(t,y)[v-r*y(1);v-r*y(2);v-r*y(3)];
options=odeset('Mass',L);
[T,Y] = ode45(fun,tspan,y0,options)

[siroos jalilyan]

Thank you; but ode45 is not true i use ode23t

[Torsten]

Either you have an old MATLAB version or your matrix L is singular.

To be on the secure side, use ode15s instead of ode45.

Best wishes

Torsten.

[siroos jalilyan]

after elapsed time (.5 second) answer trend to infinity before this time answer quite true what can i do?

[Torsten]

To give advice, you must specify the values given to the variables (L11, L12,...) .

[siroos jalilyan]

47*47 matrix

[Torsten]

Try null(L) to test whether L is singular.

And why 47x47 ? Did you modify the problem ?

[siroos jalilyan]

no because L is mutual inductance of electrical machine it have 47 DAE from beginning i try to solve 47*47 Matrix

[siroos jalilyan]

i attached answer this is a current of electrical machine plot i and try "axis([0 500 -100 100])" to show true answer

[siroos jalilyan]

if my differential equation is

v-r*i1=L11*di1/dt+L12*di2/dt+L13*di3/dt+i2*dL12/dt+i3*dL13/dt;
v-r*i2=L21*di1/dt+L22*di2/dt+L23*di3/dt+i1*dL21/dt+i3*dL23/dt;
v-r*i3=L31*di1/dt+L32*di2/dt+L33*di3/dt;i1*dL31/dt+i2*dL32/dt;

how to solve?

[Torsten]

Underdetermined system - you'll have to add six additional relations.

[siroos jalilyan]

they have not relation for additional equation actually DEA is

v-r*i1=L11*di1/dt+d(L12*i2)/dt+d(L13*di3)/dt
v-r*i2=d(L21*i1)/dt+L22*di2/dt+d(L23*i3)/dt
v-r*i3=d(L31*i1)/dt+d(L32*i2)/dt+L33*di3)/dt

above i Expansion derivatives

[Torsten]

Then the L-matrix contains given functions of time ?

[siroos jalilyan]

yes L variable with t

L=f(t)

[Torsten]

Write your system as

L11*di1/dt+L12*di2/dt+L13*di3/dt=-(i2*dL12/dt+i3*dL13/dt)+v-r*i1;
L21*di1/dt+L22*di2/dt+L23*di3/dt=-(i1*dL21/dt+i3*dL23/dt)+v-r*i2;
L31*di1/dt+L32*di2/dt+L33*di3/dt=-(i1*dL31/dt+i2*dL32/dt)+v-r*i3;

and define 'Mass' as a function instead of a constant Matrix with its entries as before.

Since the L matrix is given, you will also be able to calculate the derivatives of its elements with respect to t so that you can evaluate the right-hand side expressions.

Best wishes

Torsten.

[siroos jalilyan]

look at this code

main code

clc;
clear all;
global  DLskr nbar Lskr L LSS Lrr  Tem k T;
nbar=28;
k=0;
ChMachine;
DLskr=zeros(3,nbar);
t0=0;
tf=10;
tspan=[t0 tf];
i0=zeros([nbar+3 1]); % Current intial condtion
w0=0;                 % intial angular Velocity
ic=[i0;w0];           % Machine intitial condtion

stepsize=1e-3;
options=odeset('Mass',@Mt,'MaxStep',stepsize,'MStateDependence','none','Masssingular','yes','RelTol',10,'AbsTol',10);
[t,i] = ode15s(@wffun,tspan,ic,options);

                     %%%%%%%%%%%%%PLot %%%%%%%%%%%%%%
  figure(1) 
  plot(t,i(:,1:3))
  ylabel('Stator Current')
  xlabel('t')
  figure(2)
  plot(t,i(:,4:nbar+3))
  ylabel('Rotor Current')
  xlabel('t')
  figure(3)
  plot(t,i(:,nbar+4))
  ylabel('Angular Speed')
  xlabel('t')
  figure (4)
  plot(linspace(t0,t(length(t)),length(Tem)),Tem);
  ylabel('Tourqe')
  xlabel('t')

[siroos jalilyan]

and first function

% function ChMachin
global  L nbar  Lrr LSS Alpha delta Lm Lms p Rs Rr Vmax F J 

              %%%%%%%%Data of Macine  %%%%%%%%
  Vmax=220;                  %   Machine Voltage
  u0=4*pi*10e-7;             %   u0
  l=120e-3;                  %   Effective Length of Rotor
  g=.28e-3;                  %   Air gap Length
  N=156/3;                   %   Number of winding
  r=.012;                    %   Radiuse of Rotor
  p=2;                       %   Number of  Motors Pole 
  J=.002;                    %   Inertia of Rotor
  F=.001;                    %   Friction Coefficient                  
  Lb=0;                      %   Inductance of each Bar
  Le=0;                      %   Endig Rotor Inductance
  Alpha=2*pi/nbar;           %   Angle between two adjacent Bar 
  delta=Alpha/2;
  Rs=1.5;
  Rr=1e-2;

             %%%%%%%%Charactrist of Machine %%%%%%%%

LASAS=((N/(2*p))^2)*pi*u0*r*l/g; %  Self Inductance of Phase A Stator 
LBSBS=((N/(2*p))^2)*pi*u0*r*l/g; %  Self Inductance of Phase B Stator
LCSCS=((N/(2*p))^2)*pi*u0*r*l/g; %  Self Inductance of Phase C Stator
LASBS=-LASAS/2;                  %  Mutual Inductance between two Phase
LASCS=-LASAS/2;                  %  Mutual Inductance between two Phase
LBSCS=-LASAS/2;                  %  Mutual Inductance between two Phase
Lls=7e-3;                        %  Stator windings Leakage inductance

              %%%%%%%%Creat mutual Inductance between Stator and Rotor %%%%%%%%

Lms=(N/2)^2*pi*u0*r*l/g;
Lm=(4/(pi*N))*Lms*sin(p*delta);

for(j=1:nbar) 

      Lskr(1,j)=Lm*cos(p*((j-1)*Alpha+delta));
      Lskr(2,j)=Lm*cos(p*((j-1)*Alpha+delta)+2*pi/3);
      Lskr(3,j)=Lm*cos(p*((j-1)*Alpha+delta)-2*pi/3);

end

      Lkk=(u0*l*r/g)*Alpha*(1-Alpha)/(2*pi);    %  Self  Inductance of each Bar
      Lkn=-(u0*l*r/g)*(Alpha)^2/(2*pi);         %  Mutal Inductance of each Bar with Other

                 %%%%%% Creat Rotor Inductance Matrix  %%%%%%

   Lrr=Lkn*ones(nbar,nbar)-diag((Lkn+2*(Lb+Le))*ones(nbar,1));             %  Rotor Inductance Matrix
   Lrr=diag((Lkk+2*(Lb+Le))*ones(nbar,1))+Lrr;

                        %%%%%%%%%%%%    %%%%%%%%%%%%
  LSS=[Lls+LASAS LASBS LASCS;LASBS LASAS+Lls LASCS;LASCS LASBS Lls+LASAS]; %  Stator Inductance Matrix
  L=[LSS Lskr;Lskr' Lrr];                                                  %  Machine Inductance Matrix
  L(nbar+4,nbar+4)=J;

[siroos jalilyan]

second function

function V=wffun(t,i)
global ir  DLskr  k nbar Lm p Alpha delta Rs Rr Tem Vmax F J...
    LSS Lrr L Lskr T
k=k+1;
TL=0.0001;
a=@(t) Vmax*sin(100*t);
b=@(t) Vmax*sin(100*t-120*pi/180);
c=@(t) Vmax*sin(100*t+120*pi/180);
v=[a(t);b(t);c(t)];
%  keyboard
 thetan=(0*pi/180)*t;
for j=1:nbar 
    Lskr(1,j)=Lm*cos(p*(thetan+(j-1)*Alpha+delta));
    Lskr(2,j)=Lm*cos(p*(thetan+(j-1)*Alpha+delta)-2*pi/3);
    Lskr(3,j)=Lm*cos(p*(thetan+(j-1)*Alpha+delta)+2*pi/3); 
end
for j=1:nbar 
    DLskr(1,j)=-Lm*p*sin(p*(thetan+(j-1)*Alpha+delta));
    DLskr(2,j)=-Lm*p*sin(p*(thetan+(j-1)*Alpha+delta)-2*pi/3);
    DLskr(3,j)=-Lm*p*sin(p*(thetan+(j-1)*Alpha+delta)+2*pi/3); 
end
Tem(k)=i(1:3)'*DLskr*i(4:nbar+3);
L=[LSS Lskr;Lskr' Lrr];
L(nbar+4,nbar+4)=J;
V=[v-Rs*i(1:3)-(i(4:nbar+3)'*DLskr')';...
    (Rr*i(4:nbar+3)'-i(1:3)'*DLskr)';Tem(k)-F*i(nbar+4)];

[siroos jalilyan]

third function

function M=Mt(t)
global L 
M=L;

[siroos jalilyan]

this is analyze electrical machine code when machine is station and angular speed is 0 answer exactly true but when machine start rotating answer infinity

in function V=wffun(t,i) (second function) theatn=0; while thetan =w*t answer infinity

as you said above i do (ending line second function)

V=[v-Rs*i(1:3)-(i(4:nbar+3)'*DLskr')';...
    (Rr*i(4:nbar+3)'-i(1:3)'*DLskr)';Tem(k)-F*i(nbar+4)];

Lskr variable when machine is rotating and in any step Lskr updating

Dlskr is derivative of Lskr

|

[Torsten]

function V=wffun(t,i)
global ir  DLskr  k nbar Lm p Alpha delta Rs Rr Tem Vmax F J...
    LSS Lrr L Lskr T
k=k+1;
TL=0.0001;
a=@(t) Vmax*sin(100*t);
b=@(t) Vmax*sin(100*t-120*pi/180);
c=@(t) Vmax*sin(100*t+120*pi/180);
v=[a(t);b(t);c(t)];
%  keyboard
 thetan=(0*pi/180)*t;
 d_thetan_dt=(0*pi/180);
for j=1:nbar 
    Lskr(1,j)=Lm*cos(p*(thetan+(j-1)*Alpha+delta));
    Lskr(2,j)=Lm*cos(p*(thetan+(j-1)*Alpha+delta)-2*pi/3);
    Lskr(3,j)=Lm*cos(p*(thetan+(j-1)*Alpha+delta)+2*pi/3); 
end
for j=1:nbar 
    DLskr(1,j)=-Lm*p*d_thetan_dt*sin(p*(thetan+(j-1)*Alpha+delta));
    DLskr(2,j)=-Lm*p*d_thetan_dt*sin(p*(thetan+(j-1)*Alpha+delta)-2*pi/3);
    DLskr(3,j)=-Lm*p*d_thetan_dt*sin(p*(thetan+(j-1)*Alpha+delta)+2*pi/3); 
end
Tem(k)=i(1:3)'*DLskr*i(4:nbar+3);
L=[LSS Lskr;Lskr' Lrr];
L(nbar+4,nbar+4)=J;
V=[v-Rs*i(1:3)-(i(4:nbar+3)'*DLskr')';...
    (Rr*i(4:nbar+3)'-i(1:3)'*DLskr)';Tem(k)-F*i(nbar+4)];

[siroos jalilyan]

w=i(nbar+4)

thetan=w*t;

and try

   thetan=i(nbar+4)*t;
   if(i(nbar+4)==0)
   d_thetan_dt=1;
   else
       d_thetan_dt=i(nbar+4);
   end

[Torsten]

Your problem is that the code doesn't run ? Doesn't produce the results you expect ?

[siroos jalilyan]

while

thetan=i(nbar+4)*t;
   if(i(nbar+4)==0)
   d_thetan_dt=1;
   else
       d_thetan_dt=i(nbar+4);
   end

answer unexpect

actully thetan=0 answer is true

if motor rotate thaheta~=0

and Doesn't produce the expect results

[Torsten]

i is the vector of solution variables - thus d(i*t)/dt=t*di/dt+i.

[siroos jalilyan]

I could not understand you mean!!!!

[Torsten]

I mean that if

thetan=i(nbar+4)*t;

then

d_thetan_dt = t*d(i(nbar+4))/dt + i(nbar+4)

not what you wrote above.

[siroos jalilyan]

pay attention Tem(torque) ~= 0 i was try your code but it is not true.

in my code every thing is true until thetan=0

[siroos jalilyan]

if i dont expansion d(L*i) then

d(L11*i1)/dt+d(L12*i2)/dt+d(L13*di3)/dt = v-r*i1
d(L21*i1)/dt+d(L22*i2)/dt+d(L23*i3)/dt = v-r*i2
d(L31*i1)/dt+d(L32*i2)/dt+d(L33*i3)/dt = v-r*i3

witn condtion L11 and L22 and L33 are constant then build mass matrix and solve

how do i do?

[Torsten]

Define the system as

dy11/dt + dy12/dt + dy13/dt = v-r*i1
dy21/dt + dy22/dt + dy23/dt = v-r*i2
dy31/dt + dy32/dt + dy33/dt = v-r*i3
0 = y11-L11*i1
0 = y12-L12*i2
0 = y13-L13*i3
0 = y21-L21*i1
0 = y22-L22*i2
0 = y23-L23*i3
0 = y31-L31*i1
0 = y32-L32*i2
0 = y33-L33*i3

Twelve equations in twelve unknowns.

But I think this will result in the error message that the system has index greater than 1.