Solving a parameter estimation problem using trajectory optimization

Question

This is a follow-up to my previous question here

I've the following system of equations for studying information flow in the below graph,

$$ \frac{d \phi}{dt} = -M^TDM\phi + \text{noise effects} \hspace{1cm} (1)$$

Here, M is the incidence matrix of the graph

$\phi$ is a vector with variables [ A B C D E F].

I've solved the above odes to obtain the time series data of variables A,B,C,D,E,F.

Using the time-series data obtained from the above step, I'd like to do determine $\tilde{D}$ for the following system

$$ \frac{d \phi}{dt} = -M^T\tilde{D}M\phi \hspace{1cm} (2)$$

Note: The entries in the diagonal elements of $\tilde{D}$ are the edge weights.

In summary: Equation (1) (with noise effects) is solved using prior values of the diagonal matrix, D and the time series profiles of variables in each node are obtained. I want to determine a modified D i.e $\tilde{D}$ that can generate the same time series profile that was generated while solving equation (1).

Based on the solution provided in my previous post, I want to solve this as an optimization problem of the form $$\mathsf{K} = \int_{0}^{t_{f}} ||\phi(t) - \hat{\phi}(t)||^{2} dt$$

$$\tilde{D}, \hat{\phi}(0) = \text{argmin} \ \mathsf{K}(\tilde{D},\hat{\phi}(0)) = \text{argmin} \ \int_{0}^{t_{f}} ||\phi(t) - \exp{(-M^{T} \tilde{D} M t)} \hat{\phi}(0)||^{2} dt$$

I'd like to solve this optimization problem using fmincon in MATLAB.

The constraints will be the dynamical system presented in equation 1 above. I read through some of the procedures given in the literature and I want to use the trapezoidal rule to approximate dynamical constraints. However, I am not sure how to specify the constraints as non-linear equality constraints in MATLAB. Also, $\phi$ is a vector and I'd like to know if there is an easy way to express the constraints using the trapezoidal rule, i.e in a matrix form.

I'd also like to know if the integral form of the objective function should also be approximated using trapezoidal rule. Is it required to specify upper and lower bounds apart from the objective and equality constraints?

Any suggestions on how to proceed will be really helpful.

If there are examples for solving these kinds of problems, links to those will be useful.

EDIT: Template of implementation algorithm suggested by whpowell96

Dhat0 = %input vector 
% fun   = @objfun;
% [Dhat,fval] = fminunc(fun, Dhat0)

%% lsqnonlin
Dhat = lsqnonlin(@(Dhat) objfun(Dhat),Dhat0)


function f = objfun(Dhat)

%% Integrator settings
tspan = %tspan 
options = odeset('abstol', 1e-10, 'reltol', 1e-9);

%% generate exact solution
    phi0 = % initial condition vector
    [t, phi]  = ode15s(@(t,phi) exact(t,phi), tspan , phi0 ,options);


%% generate approximate solution

    [t, phi_tilde]  = ode15s(@(t,phi_tilde) approx(t,phi_tilde, Dhat), tspan , phi0 ,options);


%% objective function for fminunc
    % diff = (phi - phi_tilde).*(phi - phi_tilde);
    % f = sum(diff, 'all')

%% objective function for lsqnonlin
    f  = phi - phi_tilde
end

Answer 1

I am a bit confused as to your characterization of constraints. Equation $(1)$ is not a constraint. It is the model that generated the time series data you are trying to fit. You then try to find the correct parameters $\tilde{D}$ that result in equation $(2)$ matching your time series as well as possible. I would formulate the problem as the following:

1. Generate the time series data $\phi$ at some times $t_0,\dots,t_n$ using an ODE solver in MATLAB

2. Make an objective function that does the following:

   • Take in the diagonal values of $\tilde{D}$ (I believe these are the only numbers you are solving for, but that is not very clear)
   • Solve the corresponding differential equation with $\tilde{D}$ at the same time points $t_0,\dots,t_n$ to get the vector $\hat{\phi}$
   • Return the mean squared error between the $\phi$ and $\hat{\phi}$ vectors. This will approximate the integral objective functional up to $O(\Delta t)$, so it should be fine if you take enough time points.

3. Plug this new function into fminunc.

Using the ODE solver to compute $\hat{\phi}(t)$ will be much more stable than computing the matrix exponential and repeatedly multiplying. This formulation should also not take too long to run since you are only solving for 5 paramters (I think) and your ODE system is small.

Edit: lsqnonlin may be a better choice and requires a slight modification of the above advice in that you do not have to compute the mean squared error yourself. You must instead supply lsqnonlin with the vector of residuals between the two trajectories.

Answer 2

Your cost function can also be written as

$$ K = \int_0^{t_f} \left(\phi(t) - e^{-M^\top \tilde{D}\,M\,t} \hat{\phi}(0)\right)^\top \left(\phi(t) - e^{-M^\top \tilde{D}\,M\,t} \hat{\phi}(0)\right) dt. $$

When minimizing that cost function with respect to $\tilde{D}$ and $\hat{\phi}(0)$ it would be equivalent to minimizing the following cost function

$$ K = \hat{\phi}(0)^\top L_1 \hat{\phi}(0) -2\,L_2\,\hat{\phi}(0), $$

with

$$ L_1 = \int_0^{t_f} e^{-M^\top \tilde{D}^\top M\,t} e^{-M^\top \tilde{D}\,M\,t} dt, \\ L_2 = \int_0^{t_f} \phi(t)^\top e^{-M^\top \tilde{D}\,M\,t} dt. $$

Minimizing with respect to $\hat{\phi}(0)$ gives

$$ \hat{\phi}(0) = L_1^{-1} L_2^\top. $$

Substituting this back into the equivalent cost function gives

$$ K = -L_2 L_1^{-1} L_2^\top. $$

It can be noted that $L_1$ can also be obtained by solving the following Lyapunov equation

$$ M^\top \tilde{D}^\top M\,L_1 + L_1\,M^\top \tilde{D}\,M = I - e^{-M^\top \tilde{D}^\top M\,t_f} e^{-M^\top \tilde{D}\,M\,t_f}. $$

The integral of $L_2$ would still have to be evaluated. But I suspect that reducing this problem using analytical results should reduce the computation time of the cost function therefore speed up how fast this optimization problem can be solved. As already mentioned in the answer from hwpowell96 you can just use an unconstrained solver, such as fminunc or fminsearch. You do still need to provide them with a starting guess for $\tilde{D}$.