MATLAB: How to estimate or optimize the parameters of the ODE system in MATLAB 8.1 (R2013a)

Question

I have a system of ordinary differential equations (ODE) with some unknown parameters (coefficients). I want to simultaneously solve the system of differential equations as well as optimize for the unknown parameters by minimizing an objective function that depends on the solution of the ODE system.

What is the best way to do this in MATLAB?

Answer 1

There are two approaches that can be taken:

1) The optimization toolbox in MATLAB provides function such as FMINSEARCH, LSQNONLIN, FMINCON etc., which can be used for optimizing parameters while minimizing an objective function. In this case, the objective function will have to call another sub-routine which solves the differential equations using ODE solvers such as ODE23, ODE45, ODE23s, ODE113, or ODE15s. The ODE solver in turn will call the function where the differential equations are implemented.

Attached is a simple example that optimizes a single parameter for a 1 state simulation (that is a first order ODE). Note that the ‘optimization toolbox’ is needed to execute the example.

The workflow can be described as follows:

a) Choose an ODE solver (i.e. ode45 or ode15s), then write the update function for the differential equations.

b) Write an objective function that takes in the values of the parameters, solves the ODE for those particular values, and then calculates the cost function (such as the difference between the experimental and simulated data) that needs to be minimized.

c) Use an optimization function like LSQNONLIN or FMINCON to minimize the objective function.

d) Use this section of the documentation to address potential issues that go along with a numerical optimization routine wrapped around a numerical simulation:

https://www.mathworks.com/help/optim/ug/optimizing-a-simulation-or-ordinary-differential-equation.html

For example, some potential problems can occur due to the variable time-stepping used by many ODE solvers. This causes the gradient of the objective function to be ‘non-smooth’. Many optimization routines depend on computing the gradient, the change of the function value when it is perturbed. When the function value is determined by the NUMERICAL solution of an ODE, the function value may not depend smoothly on the perturbation. If you move the perturbation up a little bit, the function moves one way; if you move it up a little bit more, the function may move the opposite way. This non-smoothness is because of the variable step size used by the ODE solvers.

A good reference for this problem is Hairer, Norsett, and Wanner, Solving Ordinary Differential Equations I, Springer-Verlag, 1991, pg. 200.

One solution for the above problem is to eliminate the error checking done by an explicit ODE Solver (e.g. ode45) by forcing the solver to take fixed steps. This eliminates the need for error checking to choose a step size, but may cause accuracy problems. To do "fixed-step" simulation with the ODE Suite, you can use ODESET to turn up both the ‘RelTol’ and the ‘AbsTol’ to large values (e.g. 1), and set the ‘MaxStep’ and the ‘InitialStep’ to some fixed-step value (visit the documentation page for more information on ODESET). Note that this technique will not work with the implicit solvers such as ode15s.

To get output at specific times, you can pass in the vector of times ‘[t1 t2 ... tn]’ as the ‘tspan’ vector:

[t, y] = ode23(myfun, [t1 t2 ... tn], y0, options);

By combining these two techniques, you can avoid some of the difficulties described in Hairer, Norsett, and Wanner, but there is always the problem of not knowing the accuracy when doing fixed step simulation.

Another solution is to use PATTERNSEARCH for the optimization. It is a derivative-free solver. Note that this function belongs to the ‘Global Optimization toolbox’.

2) A second way of estimating the parameters in an ODE is to use the "Nonlinear Grey Box" models of System Identification Toolbox. This feature has been available since R2007a.

The workflow is as follows:

a) Prepare a data object containing the signal that you want to optimize your ODE coefficients for.

b) Write an M- or a MEX-file for your ODE that represents the ODE as a set of first order differential equations. This file returns the state derivatives and signal values as a function of time, forcing function (if any), coefficient values and initial conditions.

c) Then create an IDNLGREY object that encapsulates the ODE in a model form. It uses the ODE file from previous step as its main component. It also defines the initial values of coefficients, dimension of signal (in case of multi variable ODEs), names/units of signals, type of optimizer to use etc.

d) Use the estimation command PEM to optimize the ODE's coefficients to maximize the fit between the ODE's simulation results and the provided external signal. For this, a variety of optimization routines are available (line-search schemes such as Gauss-Newton, Levenburg-Marquardt, Steepest Descent and trust region approach (lsqnonlin)).

These techniques work for both linear and nonlinear ODEs. However, for linear ODEs, a special object called IDGREY is also available which facilitates a more efficient approach and provide additional functionality such as prescribing stability constraints, modeling disturbance and plotting time/frequency domain response of the model ("model" is a state-space representation of the ODE). See the demos in system Identification Toolbox for some examples.

If you're doing this type of work within a Simulink simulation, Simulink Design Optimization allows you to optimize system parameters while using a variable step size solver. It is designed to handle the numerical difficulties which arise in these situations such as the aforementioned issues when using the optimization toolbox and ode solvers directly. You can find more information on Simulink Design Optimization at:

https://www.mathworks.com/products/sl-design-optimization.html