# Solving a parameter estimation problem using trajectory optimization

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


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.

• Thanks a lot for the response. Solve the corresponding differential equation with $\tilde{D}$ at the same time points t0,…,tn to get the vector $\hat{ϕ}$ I don't think this is possible. I wan to estimate the diagonal entries of $\tilde{D}$ , these are 6 unknowns. So it's not possible to use the ode solver to estimate $\hat{ϕ}$ Or probably, I misunderstood what you are saying – Natasha Mar 18 at 2:33
• Estimation of $\tilde{D}$ is done by taking successive guesses of of $\tilde{D}$, and updating them so that the next guess lowers the error between the true and predicted trajectory. Every optimization algorithm works by evaluating a loss function as a function of parameters. It seems like you are trying to do things analytically, which may be useful, but not practical. If it is the case that $\tilde{D}$ is actually a dense matrix, then you can instead do the same optimization over all 36 components and then take the diagonal elements of the optimal solution – whpowell96 Mar 18 at 2:47
• @Natasha I think powell’s purpose is to define a function that takes a $\tilde{D}$ and solve the ODE numerically, but for actual process you need to have at least an initial guess and then iteratively follow the powell’s procedure to minimize the error or cost function. – Alone Programmer Mar 18 at 2:49
• This is what fminunc and every other algorithm have to do under the hood. Numerically solving the ODE is a faster and more stable that doing the matrix exponential formulation in the OP – whpowell96 Mar 18 at 2:50
• Unless I am drastically misunderstaing the problem, this is your best bet. I have done this with several problems and it works for systems where solving the ODEs doesn't take forever or you don't have a ton of parameters – whpowell96 Mar 18 at 2:54

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}$$.

• Thanks a lot for the response. I did try implementing what hwpowe96 suggested by providing an initial guess for $\hat{D}$ in fminunc. My first try was to solve the problem with zero noise effects in eq(1). In such a case, I expect the solver to return $\hat{D}$ = D. But this didn't work. Now, I am trying to solve this as a constrained optimization problem, providing dynamics as constraints. Please check my post here – Natasha Mar 21 at 5:15
• @Natasha I tried to implement this myself and noticed that your proposed cost function might cause issues. Namely, the contribution of the noise from time zero to current time will all influence the state at the current time. This might eventually dominate over the dynamics caused by the initial conditions. Therefore, it might be better to define the estimated state as $\hat{\phi}(t) = e^{-M^\top \tilde{D}\,M\,\Delta t} \phi(t-\Delta t)$. This also removes the need to estimate $\hat{\phi}(0)$ and hopefully reduces the effect of the noise as much as possible. – fibonatic Mar 21 at 16:16
• Thanks a lot for the response. Could you please let me know if you could obtain the initial parameter guess while setting W (noise) = 0? When W = 0 , eq(1) will be equal to eq(2) and I'd expect fminunc or lsqnonln to return $\hat{D}$ = $D$. Please note: the initial guess for $\hat{D}$ = D. – Natasha Mar 21 at 16:35
• Also, from what has been suggested by whpowell96 Using the ODE solver to compute ϕ^(t) will be much more stable than computing the matrix exponential and repeatedly multiplying. ` I tried to define a separate function for computing $\hat{\phi}$ instead of computing exponential. This function is called inside the objective function with updated parameters returned by fminunc. – Natasha Mar 21 at 16:42
• @Natasha I assume that $\phi(t)$ is only known at discrete times and if those times are equally spaced you would also have that $\Delta t$ is constant and you would only have to evaluate $e^{-M^\top \tilde{D}\,M\,\Delta t} once when evaluating the cost function. Note for these discrete times the cost function integral would be better written as a summation (otherwise you would also have to start interpolating in order to get all$\phi(t)$). But it should also work to just solve the ODE for$\Delta t$seconds with$\phi(t-\Delta t)$as initial condition in order to obtain$\hat{\phi}(t)\$. – fibonatic Mar 21 at 16:54