I am developing a practical work of the following system of ode

\begin{align}x'(t) &= k_1h(t) - (k_2+k_3)x(t)\\ y'(t) &= k_3x(t)\end{align}

and $z(t) = (1-k_4)(x(t)+y(t))+k_4h(t)$, where $h(t) = 2t+1$

We now consider the following cost function:

$$J(u) = \frac{1}{2}\int_{t_0}^{t_f}|z_u(t)-\bar{z}(t)|^2w(t)dt+\frac{\lambda}{2}|u-u_{ap}|^2.$$ Where $u=(k_1,k_2,k_3,k_4)$, the function $w(t)$ modulates the weight that is given to each value taken from the $\bar{z}$ data ($\bar{z}$ is an interpolated function of the data). If we consider that all the values have the same importance we can only take $w(t) = 1.$ The second term present in $J(u)$ is interesting when we have an idea of the values of the parameters to optimize. It therefore allows you to stay neighborhood to these values.

In order to be able to minimize $J(u)$ we introduce the Lagrangian:


where $Y(t)=(x(t),y(t))$ and $C(Y(t),u)=z_u(t)=(1-k_4)(x(t)+y(t))+k_4h(t).$

So $$J(u) = L(Y_u,u,p).$$

When $D_uJ(u) = 0$ leads to the resolution of the following adjoint equation

\begin{align}-p'(t) &= (D_yf)^*p(t)-(C(Y(t),u)-\bar{z}(t))w(t)(D_yC)\\ p(T)&=0 \end{align}

I want to minimize $J(u)$ using a constant $\rho$-step gradient method:

$$u^{k+1}=u^k-\rho D_uJ(u^k).$$

So $D_uJ(u) = D_uL(y,u,p) = \lambda(u-u_{ap})\cdot\bar{u}+\int_0^T(C(Y(t),u)-\bar{z}(t)w(t))D_uC^*-D_uf^* p(t))\cdot \bar{u}.$

Where $\bar{u}$ is a vector of $\mathbb{R}^4$ and \begin{equation*} D_yf = \begin{pmatrix} -(k_2+k_3) & 0 \\ k_3& 0 \end{pmatrix}, \end{equation*} \begin{equation*} D_uf = \begin{pmatrix} h(t) & -x(t) & -x(t) & 0 \\ 0 & 0 & x(t) & 0 \end{pmatrix}, \end{equation*} \begin{equation*} D_yC = \begin{pmatrix} 1-k_4 & 1-k_4 \end{pmatrix}, \end{equation*} \begin{equation*} D_uC = \begin{pmatrix} 0 & 0 & -(x(t)+y(t))+h(t) \end{pmatrix} \end{equation*}

I have tried in several ways to "build" the integral function but as $C(Y(u),t)$ and $\bar{z}$ are interpolated functions and in addition the matrices $D_uJ(u),\dots$ shown also contain these interpolated functions, how can I implement?

Here my code:

import numpy as np
import matplotlib.pyplot as plt
from scipy.interpolate import interp1d
from scipy.integrate import odeint

# Creation of system of ode
def myfct_R(Y, t, params):
    k1, k2, k3, k4 = params
    x, y = Y
    h = 3*t+2 # pour tester la convergence des schémas
    return np.array([k1*h-(k2+k3)*x, k3*x])
# Resolution of system of ode
def z_fct(t,y0,params):
    k1, k2, k3, k4 = params
    sol_fct = odeint(myfct_R, y0, t,args=(params,))
    return (1-k4)*(sol_fct[:,0]+sol_fct[:,1])+k4*h(t), Y

params = [0.01,0.02,0.03,0.1]
Nt = 2**14 +1
y0 = [0.0,0.0]
t = np.linspace(0, tf, Nt) # tableau de t = 0->tmax avec N valeurs
zfunct,Z = z_fct(t,y0,params)

#---------------Adjoint equation-----------------


Cyu = interp1d(t,zfunct)
def Padj_fct(p, t,  params):
    p1, p2 = p # p(t) = (p1(t),p2(t))
    k1,k2,k3,k4 = params
    w = 2*t**2 -2*t -1    
    Cyui = Cyu(t)#interp1d(t,zfunct) 

    Dyf = np.array([[-(k2+k3),0],[k3,0]]) 
    Dyfa = np.transpose(Dyf)
    DyC = np.array([1-k4,1-k4])    
    zb = h(t)    
    H = Cyui - zb
    return Dyfa.dot(p)-(H*w)*np.transpose(DyC)

def sol_adj(t,y0,params):
    k1, k2, k3, k4 = params
    Y = odeint(Padj_fct, y0, t,args=(params,))

    return Y

padj = sol_adj(t,y0,params)

And from here I don't know how to proceed with the calculation of the gradient method, I would also like to know if the Padj_fct function is good to numerically approximate the attached equation presented above.

Note: $DyC^* = DyC^T$ is the transpose matrix

  • 1
    $\begingroup$ Seems like a homework... Please don't post images instead of writing... $\endgroup$ – Alone Programmer Jan 23 at 22:49

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.