How do I estimate and optimize the parameters of my ODE system in MATLAB?

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

[MathWorks Support Team]

You can approach this problem in two ways. The first one uses the Optimization Toolbox, and the second one uses the System Identification Toolbox.

Using the Optimization Toolbox

The Optimization Toolbox in MATLAB provides optimization functions such as 'fminsearch', 'lsqnonlin', and 'fmincon'. These functions can be used for optimizing the ODE parameters by minimizing an objective function. The objective function will 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.

An example of this workflow is explained in the documentation for the Optimization Toolbox: 

https://www.mathworks.com/help/releases/R2022b/optim/ug/fit-ode-problem-based-least-squares.html

There is a small .m file attached to this article which reproduces this suggested workflow as well.

You can also consult this page to address potential issues that go along with a numerical optimization routine wrapped around a numerical simulation:

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

Some problems can appear when using variable time steps, which is a common approach by many ODE solvers. This setting causes the gradient of the objective function to be ‘non-smooth’. Many optimization routines depend on computing the gradient of the objective function, which is a measure of the change of the function value after it is perturbed slightly. 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 the book by 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. You can specify options for the ODE solution process using 'odeset'. Set both the ‘RelTol’ and the ‘AbsTol’ to large values (e.g. 1), and set the ‘MaxStep’ and the ‘InitialStep’ to a fixed-step value. Note that this technique will not work with the implicit solvers such as ode15s.

Function 'odeset' is documented here: https://www.mathworks.com/help/releases/R2022b/matlab/ref/odeset.html

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 the book by Hairer et al., but there is always the problem of not knowing the accuracy when doing fixed step simulation.

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

Using the System Identification Toolbox

A second way of estimating the parameters in an ODE is to use the "Nonlinear Grey Box" models of System Identification Toolbox. This workflow is documented in the following page:

https://www.mathworks.com/help/releases/R2022b/ident/ug/estimating-nonlinear-grey-box-models.html

The basic workflow is the following:

1. 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.
2. Create an 'idnlgrey' object that encapsulates the ODE in a model form. Use the ODE file from the previous step and provide additional information such as the initial values of coefficients, the dimension of signal (in case of multi variable ODEs), the names/units of signals, the type of optimizer to use etc.
3. Use the estimation commands 'nlgreyest' or 'pem' to optimize the ODE's coefficients. to maximize the fit between the ODE's simulation results and the provided external signal. A variety of optimization routines are available for this step, such as line-search schemes (Gauss-Newton, Levenburg-Marquardt, Steepest Descent) and trust region approach (lsqnonlin).

Cases with linear ODEs can use a special object called 'idgrey'. This object facilitates a more efficient approach and provides 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.

Additionally, if you are 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 the following page:

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