[Matlab] Problem in modeling a DAE system using ode15s

[ruwais92]

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.

 

[Greg Bernhardt]

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?

 

[ruwais92]

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!