MATLAB [Matlab] Problem in modeling a DAE system using ode15s

AI Thread Summary
The discussion centers around a coupled system of ordinary differential equations (ODEs) with an algebraic equation, specifically for a hybrid system involving fuel cells and a supercapacitor. The user shares a MATLAB function that defines the system's parameters and equations, including physical constants and equilibrium points. They express concerns about potential errors in using the ODE solver ode15s, particularly regarding the creation of a mass matrix, which is leading to diverging outputs. In response to the inquiry, another participant notes they have not pinpointed the exact error but have found a workaround that allows the system to function as intended. The exchange highlights the challenges of solving complex coupled ODEs and the importance of community support in troubleshooting coding issues.
ruwais92
Messages
2
Reaction score
0
Dear all,
I have a coupled system of ode's with an algebraic equation. The equations are as follows,
function dy=lorentz(t,y)
global Vfc;
%fcdata.m
% Data file for the hybrid system with parameters values

%***************************************
% Physical parameters
Rl = 0.024; %Value in ohm 0.07%
n = 11; % Number of fuel cells
% SuperCapacitor parameters:
Csc = 58; % Value in F Taken from owner's manual (SuperCapacitor data-sheet)
Rsc = 0.019; %Value in ohm Taken from owner's manual (SuperCapacitor data-sheet), internal DC resistance

%***************************************
% Disturbance
Iload = 25; % Value in amps
%%%%%%%%%%%%%%%%%%%%%%%%
Afc=75; % Value in cm sq
F = 96485;
mH2 = 1;

%Equilibrium point:
Tfce=353.15; % Value in K
alphafce = 0.4774; %old value = 0.4608 (old values were computed with SAILLER)
alphasce = 0.3875; %old value = 0.3875
Ifce = 35.8107; % Value in amps %old value = 36.7023
Isce = 0; % Value in amps
Vsc = 14.7; % Value in volts
Ct=17.9;

% References
Vcref = 24; % Value in volts
Lfc = 50e-6; % Value in micro H
Lsc = Lfc;
C=0.0376;
Rsc = 0.019;
T=353.15; % Value in K
Tamb=298.15; % Value in K
tau=2.06;

%%Function

x=(-1.524e-4*(y(4)))+4.110e-2;
w=(-3.185e-6*(y(4)))+4.684e-4;
z=(-1.219e-6*(y(4)))+1.108e-4;
Vfc=1.05-x.*log((y(1)*1000)/Afc)-w.*(y(1)*1000)/Afc-z.*exp(8e-3*(y(1)*1000)/Afc);
dy(1)=1/Lfc*((Vfc*n)-Rl*y(1)-(1-alphafce)*y(3));
dy(2)=1/Lsc*(Vsc-Rsc*y(2)-(1-alphasce)*y(3));
dy(3)=1/C*((1-alphafce)*y(1)+(1-alphasce)*y(2)-Iload);
dy(4)=1/Ct*(mH2*(n*y(1)/(2*F))-((Vfc*n)*((y(1)*1000)/Afc))-(Ct*((T-Tamb)/tau)));
dy=dy';
end

The function call is as follows,
% Function call for ODE %
clear all
clc
close all
global Vfc;
[t,y]=ode15s(@lorentz, [0 0.05], [35.2 0 13.4 353.15]);
subplot(231)
plot(t,y(:,1),'-');
xlabel('time');
ylabel('Ifc');
title('Ifc reponse over time');

subplot(232)
plot(t,y(:,2),'-');
xlabel('time');
ylabel('Isc');
title('Isc reponse over time');

subplot(233)
plot(t,y(:,3), '-');
xlabel('time');
ylabel('Vc');
title('Vc reponse over time');

subplot(234)
plot(t,y(:,4),'-');
xlabel('time');
ylabel('T');
title('T reponse over time');

subplot(235)
plot(t,Vfc,'-');
xlabel('time');
ylabel('Vfc');
title('Vfc reponse over time');

I believe there is an error in the usage of ode15s. I don't know how to create the mass matrix, what it means or how I'm supposed to use it. I am getting diverging outputs. Could anybody please help me?
Thank you.
 
Physics news on Phys.org
I'm sorry you are not generating any responses at the moment. Is there any additional information you can share with us? Any new findings?
 
Dear Greg,
I have not yet been able to identify the exact nature of the error in that code. However, I did manage to find a workaround. It's working fine just like the original was intended to.
Thankyou very much for your your concern. Much appreciated!
 

Similar threads

Replies
8
Views
2K
Replies
4
Views
1K
Replies
5
Views
2K
Replies
4
Views
2K
Replies
9
Views
5K
Replies
4
Views
3K
Replies
1
Views
2K
Replies
4
Views
2K
Back
Top