# How to minimize a integral function using a constant step gradient method in Python?

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:

$$L(y,u,p)=\frac{1}{2}\int_{t_0}^{t_f}|C(Y(t),u)-\bar{z}(t)|^2w(t)dt+\frac{\lambda}{2}|u-u_{ap}|^2+\int_{t_0}^{t_f}(Y'(t)-f(Y(t),u))p(t)dt,$$

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)

#------------****************------------

Cyu = interp1d(t,zfunct)
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)


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