Turn on suggestions

Auto-suggest helps you quickly narrow down your search results by suggesting possible matches as you type.

Showing results for

- Home
- /
- Analytics
- /
- Stat Procs
- /
- Re: parameter estimation in systems of ordinary differential equations

Options

- RSS Feed
- Mark Topic as New
- Mark Topic as Read
- Float this Topic for Current User
- Bookmark
- Subscribe
- Mute
- Printer Friendly Page

- Mark as New
- Bookmark
- Subscribe
- Mute
- RSS Feed
- Permalink
- Report Inappropriate Content

Posted 02-27-2022 07:27 PM
(665 views)

```
Hey guys. I'm trying to estimate parameters in a system of ordinary differential equations.
Here is the code I wrote (I know, is a disaster), with the error messages. Can anyone tell me how to fix the code?
```

title 'ODE BATALIM';

data BATALIM;

input hour V CELL S P;

datalines;

0.0 1.5 61.64 0. 30.39

1.0 2.2 34. 43. 31.18

2.0 2.9 36.34 56.06 35.95

3.0 3.6 21.92 66.29 38.87

4.0 4.3 17.75 69.12 42.31

5.0 5.0 17.01 70.14 44.17

;

proc model data=BATALIM;

dependent V CELL S P;

parms Fa mimax Ks Ki Pmax CELLmax n m;

/*Next I try to provide the initial conditions of the differential equations*/

if (hour=0) then

do;

V=1.5; CELL=61.64; S=0; P=30.39;

end;

else

do;

/*The system of differential equations is described below*/

dert.V = Fa;

dert.CELL = (mimax*S/(Ks+S+S**2/Ki)*(1-P/Pmax)**n*(1-CELL/CELLXmax)**m)*CELL - Fa*CELL/V;

dert.S = Fa*(213-S)/V - (mimax*S/(Ks+S+S**2/Ki)*(1-P/Pmax)**n*(1-CELL/CELLmax)**m)*CELL/0.017;

dert.P = (mimax*S/(Ks+S+S**2/Ki)*(1-P/Pmax)**n*(1-CELL/CELLmax)**m)*CELL*24.94 - Fa*P/V;

end;

/*

`Then I try to insert the initial values of the parameters and get the fitted equations...`

*/

fit V start=(Fa 0.5) / time=hour;

fit CELL start=(mimax 0.01 Ks=2 Ki=5 Pmax=40 n=2 CELLmax= 60 m=2) / time=hour;

fit S start=(Fa 0.5 mimax 0.01 Ks=2 Ki=5 Pmax=40 n=2 CELLmax= 60 m=2)/ time=hour;

fit P start=(mimax 0.01 Ks=2 Ki=5 Pmax=40 n=2 CELLmax= 60 m=2)/ time=hour;

run;

Error messages

WARNING: Missing value encountered for initial dert.CELL at time = 0. This variable is set to zero.

WARNING: Missing value encountered for initial dert.S at time = 0. This variable is set to zero.

WARNING: Missing value encountered for initial dert.P at time = 0. This variable is set to zero.

NOTE: At OLS Iteration 0 CONVERGE=0.001 Criteria Met.

WARNING: Missing value encountered for initial dert.CELL at time = 0. This variable is set to zero.

WARNING: Missing value encountered for initial dert.S at time = 0. This variable is set to zero.

WARNING: Missing value encountered for initial dert.P at time = 0. This variable is set to zero.

WARNING: Missing value encountered for initial dert.CELL at time = 0. This variable is set to zero.

WARNING: Missing value encountered for initial dert.S at time = 0. This variable is set to zero.

WARNING: Missing value encountered for initial dert.P at time = 0. This variable is set to zero.

And more...

9 REPLIES 9

- Mark as New
- Bookmark
- Subscribe
- Mute
- RSS Feed
- Permalink
- Report Inappropriate Content

I am probably answering the wrong question here, but I noticed that the derivative of CELL at hour=0 is, in turn, equal to zero. That may be leading to the missing values. Looking at the example on fitting systems of ODEs in the PROC MODEL documentation, I see that it is set up to only calculate derivatives if 'time ne 0' and uses a quadratic solution for time=0. I think that particular approach is designed to avoid the warnings that you are getting with your code.

The other issue is what you copied from the log - it appears that the process is "stuck" at hour=0, but I can't tell. The repetitive nature makes me think that the same error is generated at each iteration.

Given these you ,might want to reparameterize your system of equations, if at all possible. The derivatives look unwieldy, so I suspect the partials with respect to the parameters may be even worse.

SteveDenham

- Mark as New
- Bookmark
- Subscribe
- Mute
- RSS Feed
- Permalink
- Report Inappropriate Content

**Hello Steve, how are you? Thanks for responding to my help request.**

"I am probably answering the wrong question here, but I noticed that the derivative of CELL at hour=0 is, in turn, equal to zero. That may be leading to the missing values... I see that it is set up to only calculate derivatives if 'time ne 0' ..."

**This is a problem. In various physical phenomena dy/dt (x = 0) =0.**

"The other issue is what you copied from the log - it appears that the process is "stuck" at hour=0, but I can't tell. The repetitive nature makes me think that the same error is generated at each iteration."

**I agree with you. The problem is that I can't fix it.**

"Given these you ,might want to reparameterize your system of equations, if at all possible. The derivatives look unwieldy, so I suspect the partials with respect to the parameters may be even worse."

**I was also suspicious of the equations. However, I have consulted the literature and they are correct. I've tried everything here and I don't know how to solve it.**

**I will try to find a SAS code for solving EDOS systems that works and test it with my data.**

- Mark as New
- Bookmark
- Subscribe
- Mute
- RSS Feed
- Permalink
- Report Inappropriate Content

*> I was also suspicious of the equations. However, I have consulted the literature and they are correct.*

Great. Please provide the reference to the literature that you consulted. The paper, article, or blog post will tell us what problem you are trying to solve, and we might be able to find the source of your problem.

- Mark as New
- Bookmark
- Subscribe
- Mute
- RSS Feed
- Permalink
- Report Inappropriate Content

Hi Rick!

The literature I am consulting is in Portuguese. In addition, it is very long, as it is a master's thesis. However, I have an article (unfortunately, also in Portuguese) with only 15 pages that portrays exactly what I want to do: simulate fermentation processes, both in batch and fed batch. If you want, I can send it. I'll understand if you don't want to.

- Mark as New
- Bookmark
- Subscribe
- Mute
- RSS Feed
- Permalink
- Report Inappropriate Content

Unfortunately, I would not be able to understand it. I wish you good luck in solving your problem.

- Mark as New
- Bookmark
- Subscribe
- Mute
- RSS Feed
- Permalink
- Report Inappropriate Content

I am not an expert on PROC MODEL, but I noticed two things:

1. You use the parameter CELLXmax, which is not defined. Presumably, that is a typo and you meant CellMax.

2. The parameter values for CellMax and PMax are less than the data values. From the names, I assume that Cell < cellMax and P < PMax for all times? If so, change those initial parameter values:

```
parm Fa 0.5 mimax 0.01 Ks 2 Ki 5
Pmax 50 CELLmax 65 n 2 m 2;
```

3. Do you need to add a BOUNDS statement? It seems like a lot of the parameters converge to large negative numbers. Check the model to see if it constrains any parameters.

4. I don't think you need the START= options on the FIT statements, but I don't really know what you are trying to accomplish.

Sometimes getting something that runs is the first step to debugging other problems. Here is some code that runs, but it does not provide a good fit to the data. I encourage you to check the ODE formulas (are all the signs correct?).

```
title 'ODE BATALIM';
data BATALIM;
input hour V CELL S P;
datalines;
0.0 1.5 61.64 0. 30.39
1.0 2.2 34. 43. 31.18
2.0 2.9 36.34 56.06 35.95
3.0 3.6 21.92 66.29 38.87
4.0 4.3 17.75 69.12 42.31
5.0 5.0 17.01 70.14 44.17
;
proc model data=BATALIM plots=none;
dependent V 1.5 CELL 61.64 S 0. P 30.39;
parm Fa 0.5 mimax 0.01 Ks 2 Ki 5
Pmax 50 CELLmax 65 n 2 m 2;
/* do you need to add BOUNDS????? */
Bounds 45 < PMax , 65 < Cellmax;
if (hour=0) then
do;
V=1.5; CELL=61.64; S=0; P=30.39;
end;
else
do;
/*The system of differential equations is described below*/
dert.V = Fa;
dert.CELL = (mimax*S/(Ks+S+S**2/Ki)*(1-P/Pmax)**n*(1-CELL/CELLmax)**m)*CELL - Fa*CELL/V;
dert.S = Fa*(213-S)/V - (mimax*S/(Ks+S+S**2/Ki)*(1-P/Pmax)**n*(1-CELL/CELLmax)**m)*CELL/0.017;
dert.P = (mimax*S/(Ks+S+S**2/Ki)*(1-P/Pmax)**n*(1-CELL/CELLmax)**m)*CELL*24.94 - Fa*P/V;
end;
fit V CELL S P / time=hour out=ModelOut;
run;
/* visualize the results */
data results;
merge BATALIM(rename=(V=VObs CELL=CellObs S=SObs P=Pobs)) ModelOut;
by hour;
run;
ods graphics/width=400px height=300px;
proc sgplot data=results;
scatter x=Hour y=VObs;
series x=Hour y=V;
run;
proc sgplot data=results;
scatter x=Hour y=CellObs;
series x=Hour y=Cell;
run;
proc sgplot data=results;
scatter x=Hour y=SObs;
series x=Hour y=S;
run;
proc sgplot data=results;
scatter x=Hour y=PObs;
series x=Hour y=P;
run;
```

- Mark as New
- Bookmark
- Subscribe
- Mute
- RSS Feed
- Permalink
- Report Inappropriate Content

Hello how are you?

First of all, thank you very much for your help.

I tried to write a code to estimate the parameters of a system of ordinary differential equations. The result was that mess you saw.

The parameter values

Fa 0.5 mimax 0.01 Ks 2 Ki 5 Pmax 50 CELLmax 65 n 2 m 2

they are only initial estimates, since the estimation method is iterative. Pmax is expected to be close to the largest value of P in the dataset (slightly larger or slightly smaller). The same goes for CELLmax.

I used a system of equations taken from a master's thesis. Perhaps the author made a typo. I will check them by consulting the original source.

Thanks for writing and sending me code that works!

- Mark as New
- Bookmark
- Subscribe
- Mute
- RSS Feed
- Permalink
- Report Inappropriate Content

I'd like to add a bit more perspective to the comments from @Rick_SAS that I think might assist as you debug this program.

If you take the program as supplied by @Rick_SAS but make the following modifications:

- Remove the BOUNDS statement
- Set the value of PMAX to 40 (as in the original code)
- Set the value of CELLMAX to 60 (as in the original code)

Then you will receive exactly the same error ("Missing value encountered for intial dert.CELL at time = 0..." And so on.) I'm making an educated guess here (I do not know much about PROC MODEL) but I am quite certain this is the reason:

As @Rick_SAS notes, the values of CELL and P in the data exceed the values of the parameters CELLMAX and PMAX. In the most general case, this need not necessarily be a problem. (For example, a measured value might exceed a theoretical maximum because of random error during the measurement process.) BUT in your case this is almost certainly a problem.

To see why, notice the structure of the derivative equations (for example, dert.CELL). The equation contains the term "(1 - P/Pmax)**n". If P>Pmax, then the term inside the parentheses is negative. This is OK when n=2 exactly (as you specify for the initial value of the parameter). But this is not defined when, say, n=2.1. This latter fact is the reason for the error. For a nonlinear least squares problem, the derivatives of the objective function with respect to the parameters (these are computed internally by PROC MODEL) must be well-defined in an open region. For your case, the derivatives are defined only at a point (n=0).

Depending on the application, you might consider fixing the values of the parameters m and n.

- Mark as New
- Bookmark
- Subscribe
- Mute
- RSS Feed
- Permalink
- Report Inappropriate Content

`Hello. I found your remarks interesting. `

However, I did what you suggested, but fitting the curves to the experimental data remained unacceptable...

**Available on demand!**

Missed SAS Innovate Las Vegas? Watch all the action for free! View the keynotes, general sessions and 22 breakouts on demand.

What is ANOVA?

ANOVA, or Analysis Of Variance, is used to compare the averages or means of two or more populations to better understand how they differ. Watch this tutorial for more.

Find more tutorials on the SAS Users YouTube channel.