Solve DAEs using IDASolve function of SUNDIALS 2.6.2 version

조회 수: 3 (최근 30일)
Ajinkya
Ajinkya 2025년 4월 22일
댓글: Ajinkya 2025년 5월 1일
I am trying to solve the system of differential algebraic equations which are used to model the IEEE 9 bus system. The differential equations determines the dynamic behaviour of the generators while the algebraic equations consists of the stator and the network (power flow) equations.
I have attached a MATLAB code in which the function "fun_DAE" gives out all of the equations in the form of a column vectors. ' F_D ' gives the 21 differential equations (7 for each generator), ' F_SA ' gives 18 network equations (real and reactive power balances for 9 buses) and ' F_SA' gives the 6 stator equations (2 for each generator). ' F_SA' and ' F_N ' constitutes the algebraic constraints.
the differential variables are for each generator. The initial values of this state variables are calculated and stored as matrix.
I need some help in using the IDAInit, IDASetOptions and IDASolve functions or any other relevant functions to be used to get a solution.
Thank You.
  댓글 수: 34
이전 댓글 32개 표시
Torsten
Torsten 2025년 5월 1일
Implementing such a complicated DAE system does not mean writing the equations, clicking on the RUN button and getting results as expected. It's an iterative process: you will have to check the results, make a guess what could be wrong in the code when you look at the curves, check and perhaps modify parameters and equations again, run the code, ... .
I did what I could to make the code work technically. Now it's up to you to make it give results that physically make sense.
Ajinkya
Ajinkya 2025년 5월 1일
Surely. I will do the needful.
Thank you so much for your valuable suggestions and help.
I was stuck with this problem for days and now I have done something fruitful.
Grateful for the time you have given !!

댓글을 달려면 로그인하십시오.

이 질문에 답변하려면 로그인하십시오.

채택된 답변

Steven Lord
Steven Lord 2025년 4월 22일
Since you're using release R2024b, you have access to three of the SUNDIALS solvers via the ode object. This functionality was added in release R2024a as stated in the Release Notes here. Try using the object to set up your problem and solve it using "idas" as the value of the Solver property.
  댓글 수: 2
Ajinkya
Ajinkya 2025년 4월 23일
could you please guide me how do I proceed in this?
Torsten
Torsten 2025년 4월 23일
Ran in:
F = ode;
F.ODEFcn = @(t,y) 2*t;
F.InitialValue = 0;
F.Solver = 'idas';
F
F =
ode with properties: Problem definition ODEFcn: @(t,y)2*t InitialTime: 0 InitialValue: 0 Jacobian: [] EquationType: standard Solver properties AbsoluteTolerance: 1.0000e-06 RelativeTolerance: 1.0000e-03 Solver: idas Show all properties
sol = solve(F,0,10);
plot(sol.Time,sol.Solution)
●

댓글을 달려면 로그인하십시오.

추가 답변 (1개)

Ajinkya
Ajinkya 2025년 5월 1일
clear
clc
format short
format compact
% LOADS MAT FILE CONTAINING MACHINE DATA "dataM" and EXCITER DATA "dataE"
load data
% THE ABOVE DATA IS DEFINED IN THE FILE "Machine_Exciter_Data.m"
% YBUS AND LOAD FLOW PROGRAM
[Y,Vb] = YBus();
I = Y*Vb;
V = abs(Vb); th = angle(Vb);
% POWER INJECTIONS AT BUSES
Sbus = Vb.*conj(I);
% CREATING A STRUCTURE NAMED 'BUS' TO STORE VALUES OF V, th, Y
save Bus V Y th
%%
% MACHINE DATA
[D_M, H, M, Xd, Xdd] = deal(data.D_M, data.H, data.M, data.Xd,data.Xdd);
[Xq, Xqq, Td0, Tq0] = deal(data.Xq, data.Xqq, data.Td0, data.Tq0);
D = D_M.*M;
ws = 2*pi*60;
% EXCITER DATA
[KA, TA, KE, TE, KF, TF] = deal(data.KA, data.TA, data.KE, data.TE, data.KF, data.TF);
% DETERMINING THE INITIAL VALUES OF THE STATE VARIABLES
for k = 1:3
Pg(k) = real(Sbus(k));
Qg(k) = imag(Sbus(k));
Ig(k) = conj((Pg(k) + 1i*Qg(k))/(V(k)*exp(1i*th(k))));
delt(k) = angle(V(k)*exp(1i*th(k)) + 1i*Xq(k)*Ig(k));
Idq(k) = abs(Ig(k))*exp(1i*(angle(Ig(k))-delt(k) + pi/2));
Id(k) = real(Idq(k)); Iq(k) = imag(Idq(k));
Vdq(k) = V(k)*exp(1i*(th(k) - delt(k) + pi/2));
Ed(k) = real(Vdq(k) - Xqq(k)*imag(Idq(k)));
Eq(k) = imag(Vdq(k)) + Xdd(k)*real(Idq(k));
Efd(k) = imag(Vdq(k)) + Xdd(k)*real(Idq(k)) + (Xd(k) - Xdd(k))*real(Idq(k));
Rf(k) = KF(k)*(imag(Vdq(k)) + Xdd(k)*real(Idq(k)) + (Xd(k) - Xdd(k))*real(Idq(k)))/TF(k);
Vr(k) = (KE(k) + 0.0039*exp(1.555*Efd(k)))*Efd(k);
Vref(k) = V(k) + Vr(k)/KA(k);
Tm(k) = Ed(k)*real(Idq(k)) + Eq(k)*imag(Idq(k)) + (Xqq(k) - Xdd(k))*real(Idq(k))*imag(Idq(k));
w(k) = ws;
end
%%
for i = 1:3 % STATE VARIABLE VECTOR
% 1 2 3 4 5 6 7
x(:,i) = [Eq(i); Ed(i); delt(i); w(i); Efd(i); Rf(i); Vr(i)];
end
% InitVal STRUCT TO STORE THE INITIAL VALUES OF THE DIFFERENTIAL VARIABLE
y0 = [x(:,1);x(:,2);x(:,3);Id(1,:)';Iq(1,:)';V(:,1);th(:,1)];
Init.x = x;
Init.Idq = [Id;Iq];
Init.Vref = Vref; Init.Tm = Tm;
Init.V = V;
Init.th = th;
Init.Y = Y;
save y y0
save Init Init % SAVING THE INITIAL VALUE STRUCT IN DATA FORMAT
% Eq1 = y(1); Eq2 = y(8); Eq3 = y(15);
% Ed1 = y(2); Ed2 = y(9); Ed3 = y(16);
% delta1 = y(3); delta2 = y(10); delta3 = y(17);
% omega1 = y(4); omega2 = y(11); omega3 = y(18);
% Efd1 = y(5); Efd2 = y(12); Efd3 = y(19);
% Rf1 = y(6); Rf2 = y(13); Rf3 = y(20);
% Vr1 = y(7); Vr2 = y(14); Vr3 = y(21);
%
% Id1 = y(22); V1 = y(28); theta1 = y(37);
% Id2 = y(23); V2 = y(29); theta2 = y(38);
% Id3 = y(24); V3 = y(30); theta3 = y(39);
% Iq1 = y(25); V4 = y(31); theta4 = y(40);
% Iq2 = y(26); V5 = y(32); theta5 = y(41);
% Iq3 = y(27); V6 = y(33); theta6 = y(42);
% V7 = y(34); theta7 = y(43);
% V8 = y(35); theta8 = y(44);
% V9 = y(36); theta9 = y(45);
%% SOLVING THE DAE
F = ode;
F.ODEFcn = @(t,y) func(t,y,data,Init);
F.InitialValue = y0;
F.MassMatrix = massMatrix();
F.RelativeTolerance = 1e-03;
F.Solver = 'ode15s';
%F.Jacobian = @(t,y,data) jcb;
% F.SolverOptions.MaxOrder = 3; % NOT FOR ODE23x
F.SolverOptions.MaxStep = 0.1;
S1 = solve(F,0,100)
% writematrix(S1.Time,'result.xlsx','Sheet','Time');
% writematrix(S1.Solution,'result.xlsx','Sheet','State')
%%
plot(S1.Time, S1.Solution(1, :), '-r', 'DisplayName', 'Eq1');
hold on;
plot(S1.Time, S1.Solution(8, :), '-g', 'DisplayName', 'Eq2');
plot(S1.Time, S1.Solution(15, :), '-b', 'DisplayName', 'Eq3');
hold off;
xlabel('Time');
ylabel('Flux Linkages');
legend show;
grid on;

카테고리

Help Center 및 File Exchange에서 Particle & Nuclear Physics에 대해 자세히 알아보기

태그

제품


릴리스

R2024b

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by