Unwanted negative solution with ode solver

6 views (last 30 days)
Harry Parker
Harry Parker on 24 Feb 2021
Commented: Harry Parker on 26 Feb 2021
Hi,
Any help would be very much appreciated. I'm copying a report which describes the change of microlagae under multiple variables. Unfortunatley some of the outputs I'm recieving are unfeasible. Specifically the Starch(S), Lipid(L), and Active Biomass(x) solutions. I recieve negative values for the starch and Lipid which of course is not feasabile and my active Biomass is greater then my total biomass(X), again not possible. In the code it can be shown the the active biomass is just the Biomass minus the starch and Lipids( hence dxdt(2) =dxdt(1)-dxdt(3) -dxdt(4) ). All the other results produced ( for the Biomass, nitrogen, Nitrogen quota, and Acetate ) are the same as in the model in the paper, therefore none of the equations are variables have been inputted incorrectly. I've been stuck on this for a couple days and any suggestions would be greatly appreciated.
function miroalgal_modelfgt()
close all
clear
clc
%ODE function
Ao = 0.42; %Initial acetate concentration
No = 0.3824; %Initial Nitrogen concentration
qno = 0.876; %Minimum nitrogen quota gNgC-1
tspan =[0,200]; % Experiment length hours
X0 = [0.001, 0 , 0, 0, Ao, 7, No, qno]; %[X, x*, S, L, A , H, N, qn]
[t,X] = ode15s(@microchange, tspan, X0);
save('t.mat','t');
save('X.mat','X')
%Plots of given variables vs time
subplot(2,4,1);
plot(t,X(:,1));
title('Biomass');
xlabel('Time ( h )');
ylabel('Biomass ( g / L )');
subplot(2,4,7);
plot(t,X(:,2), 'r');
title('Active Biomass');
xlabel('Time ( h )');
ylabel('Active Biomass ( g / L )');
subplot(2,4,5);
plot(t,X(:,3),'r');
title('Starch');
xlabel('Time ( h )')
ylabel('Starch ( g / L )');
subplot(2,4,6);
plot(t,X(:,4), 'r');
title('Lipids');
xlabel('Time ( h )');
ylabel('Lipid ( g / L )');
subplot(2,4,4);
plot(t,X(:,5));
title('Acetate');
xlabel('Time ( h )');
ylabel('Aetate ( g / L )');
subplot(2,4,8);
plot(t,X(:,6), 'r');
title('pH');
xlabel('Time ( h )');
ylabel('pH ( g / L )');
subplot(2,4,2);
plot(t,X(:,7));
title('Nitrogen');
xlabel('Time ( h )');
ylabel('Nitrogen ( g / L )');
subplot(2,4,3);
plot(t,X(:,8));
title('Nitrogen Quota');
xlabel('Time ( h )');
ylabel('Nitrogen Quota ( g / L )');
%Function for dynamic variables
function dxdt = microchange( ~ , var)
X = var(1); %Biomass
x = var(2); %Active biomass
S = var(3); %Starch
L = var(4); %Lipid
A = var(5); %Acetate
N = var(7); %Nitrogen
qn = var(8); %nitrogen quota
%Parameters
umax = 0.106; %maximum specific growth rate h-1
ksa = 1.789; %Acetate saturation constant gCL-1
kia = 0.109; %Acetate inhibition constant gCL-1
ksi = 1.4; %Light saturation constant umolm-2s-1
kii = 186.52; %Light inhibition constant umolm-2s-1
yxa = 0.059; %Aceatate yield coeffcient gCgC-1
o = 1; %Light attenuation coefficient LgC-1m-1
pnmax = 40.445; %Maximum N uptake rate gNfC-1h-1
k = 0.3125; %Nitrogen saturation constant gNL-1
n = 18.183; %Shape parameter
on = 137.455; %Uptake regulation coeffcient LgC-1
ksn = 0.162; %Nitrogen saturation uptake coeffcient gNL-1
kin = 0.113; %Nitrogen uptake inhibition constant gNL-1
ksan = 1.004; %Acetate saturation uptake coeffcient gNL-1
kian = 1.098; %Acetate uptake inhibition constant gCL-1
r1 = 0.0420; %Rate of reaction for R1 gCgC-1
r2 = 0.1620; %Rate of reaction for R2 gNgCh-1
r3 = 0.0041; %Rate of reaction for R3 gCgC-1
r4 = 0.0049; %Rate of reaction for R4 gNgC-1h-1
kis = 0.2079; %Inhibiton constant for R1 gNL-1
ns = 3.6205; %Shape parameter for R1
k1 = 0.1771; %Regulation constant for R1
os = 0.6675; %Regulation coeffcient for R1 LgC-1
ksl = 0.0227; %Saturation constant for R3 gNL-1
kil = 0.0861; %Inhibition constant for R3 gNL-1
nl = 1.8117; %Shape parameter for R3
k2 = 0.2135; %Regulation constant for R3
kh = 4.653; %pH coeffcient LgC-1h-1
io = 125; %Incident light intensity
z = 0.25; % Distance to culture from light
I = io * exp(-o*X*z); %Beer Lambert Correlation
wh = (A / ksa)/((A / ksa) + (I / ksi)); %heterotrophic weigthing function
wi = (I / ksi)/((A / ksa) + (I / ksi)); %Phototrophic Weighting function
uha = A / (A + ksa + (A^2/kia)); %Heterotrophic rate
uii = I / (I + ksi + (I^2/kii)); %Phototrophic rate
ummax = umax * (wh * uha + wi *uii); %Maximum specific growth rate under mixotrohic conditions
u = ummax * (1 - (qno/qn))+eps; %Specific growth rate
%Initial nitrogen concentration
pnnmax = pnmax * No^n/(No^n + k^n)*exp(-on*X); %Maximum nitrgoen uptake rate
pn = pnnmax * (N / (N +ksn +(N^2/kin))) * (A / (A + ksan + (A^2 / kian))); %Nitrogen uptake rate
ni = qn * X; %Internal nitrogen concentration
Aint = Ao - A; %Internal carbon concentration
R1 = r1 * (ni^ns / (ni^ns +(ni^2/kis)^ns))*(k1 / (k1 + (N/No))) * (1 + (1/u) * exp(os*Aint)) * u * x; %Starch synthesis rates
R3 = r3 * (ni^nl / (ni^nl + ksl^nl + (ni^2/kil)^nl))*(k2 / (k2 + (N/No))) * (1 + u) * x; % Lipids Synthesis rates
R2 = (r2/qn)*X; %Starch degradation
R4 = (r4/qn)*X; %Lipid segradation
dxdt = zeros(8,1);
%dXdt
dxdt(1) = u * X; %Biomass
%dxdt
dxdt(2) = u * X - R1 + R2 - R3 + R4; %Active Biomass
%dSdt
dxdt(3) = R1 - R2; %Starch
%dLdt
dxdt(4) = R3 - R4; %Lipid
%dAdt
dxdt(5) = -(1/yxa)*(uha/(uha+uii))*dxdt(1); %Acetate
%dhdt
dxdt(6) = kh * dxdt(2); % pH
%dndt
dxdt(7) = -pn * X; %Nitrogen
%dqndt
dxdt(8) = pn - u * qn; %Nitrogen quota
end
end

Answers (1)

Steven Lord
Steven Lord on 24 Feb 2021
Either:
  • specify the NonNegative option using odeset to control which components of the solution are not allowed to become negative or
  • use an events function to stop the solver when one of the components would become negative, make the appropriate changes to the function and/or the value of the solution, and use the solution from the time when the solver stopped (modified if necessary) as the initial condition for another call to the solver starting at the time the first call stopped.
See the ballode example for an illustration of the latter technique.

Categories

Find more on Large Files and Big Data in Help Center and File Exchange

Community Treasure Hunt

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

Start Hunting!