Optimizing the ratio between two differential equations

Question

I have two second-order differential equations describing two systems, both with the same form but different coefficients:

$$ x_v(t) = A_v \frac{d^2u(t)}{dt^2} + B_v \frac{du(t)}{dt} + C_v u(t) - D_v \frac{dx_v(t)}{dt} - \frac{d^2 x_v(t)}{dt^2} $$

$$ x_h(t) = A_h \frac{d^2u(t)}{dt^2} + B_h \frac{du(t)}{dt} + C_h u(t) - D_h \frac{dx_h(t)}{dt} - \frac{d^2 x_h(t)}{dt^2} $$

And I'm trying to select the control function $u(t)$ to maximize one system while minimizing the other; maximize $\int{x_v\ dt}$ (or alternatively the energy integral) while making $x_h$ as close to 0 as possible.

I've tried various optimal control and MPC techniques like iLQR (just makes a step function), pseudospectral optimizers like PSOPT (produces some sensible result with an extremely fine node resolution, but extremely computationally demanding and I can do better by guessing), but I feel like there must be something I'm missing - something about the ratio or the DE that would simplify the problem...

If one system is simply set to a reasonable output and the other to 0, an overdetermined system is obtained, which does not seem super useful.

This problem seems to be similar to a class called "simultaneous stabilization" or "underactuated systems", but the methods don't seem to be particularly applicable.

If anyone has any examples of an equivalent problem or a solution technique, I would be grateful! Thanks!

Ed. per questions:

All coefficients can be assumed to be positive. Both sets of coefficients are similar (whatever that's worth). The time horizon of interest is positive and finite.

The physical context of the problem is called cellular electroporation: this is the paper with transfer functions I'm working off of.

Answer 1

As

$$ \cases{ X_v(s) = \frac{A_v s^2+B_v s+ C_v}{s^2+sD_v + 1}U(s) + \frac{\dot x_v(0)+ s x_v(0)}{s^2+sD_v + 1}\\ X_h(s) = \frac{A_h s^2+B_h s+ C_h}{s^2+sD_h + 1}U(s) + \frac{\dot x_h(0)+ s x_h(0)}{s^2+sD_h + 1} } $$

we have

$$ \cases{ x_v(t) = \int \Phi_v(\tau)u(t-\tau)d\tau+\mathcal{L}^{-1}\left[\frac{\dot x_v(0)+ s x_v(0)}{s^2+sD_v + 1}\right]\\ x_h(t) = \int \Phi_h(\tau)u(t-\tau)d\tau+\mathcal{L}^{-1}\left[\frac{\dot x_h(0)+ s x_h(0)}{s^2+sD_h + 1}\right]\\ } $$

Here $\Phi_v(\tau), \Phi_h(\tau), \mathcal{L}^{-1}\left[\frac{\dot x_v(0)+ s x_v(0)}{s^2+sD_v + 1}\right], \mathcal{L}^{-1}\left[\frac{\dot x_h(0)+ s x_h(0)}{s^2+sD_h + 1}\right]$ can be obtained in closed form. The calculations for $\int \Phi_v(\tau)u(t-\tau)d\tau, \int \Phi_h(\tau)u(t-\tau)d\tau$ can also be obtained in closed form with the help of a symbolic processor, after proposing $u(t)$ as a cubic spline over a grid $\{t_0, t_1,\cdots, t_n\}$ of base points, being then $u(t)= u(u_0, u_1,\cdots, u_n,t)$.

Finally we can solve the maximization problem

$$ \max_{u_0,u_1,\cdots,u_n}\int_T (x_v^2(t,u_0,u_1,\cdots,u_n)-\rho^2 x_h^2(t,u_0,u_1,\cdots,u_n)) dt\ \ \text{s. t.}\ \ \|u\| \le u_{max} $$

Follows a simplified MATHEMATICA script assuming $A_v=3,B_v=1,A_h=0.1,B_h =1,C_v=1,C_h=1,D_v=1,D_h=2, x_v(0) = \dot x_v(0)=x_h(0) = \dot x_h(0)=1$ and $u$ as a linear spline.

ClearAll["Global`*"]
parms = Thread[{Dv, Dh, Av, Ah, Bv, Bh, Cv, Ch, uv0, uh0, duv0, 
 duh0} -> {1, 2, 1/2, 3, 1/10, 1, 1, 1, 1, 1, 1, 1}];
tmin = 0;
tmax = 3;
n = 5;
tau0 = 100;
dt = (tmax - tmin)/n;
For[k = 1, k <= n + 1, k++, t0[k] = tmin + dt (k - 1)];
nt = n + 1;
U0 = Table[u0[k], {k, 1, nt}];
ut[u1_, u2_, t1_, t2_, t_] := u1 + (u2 - u1)/(t2 - t1) (t - t1)
uf[t_] := Piecewise[
Table[{ut[u0[k], u0[k + 1], t0[k], t0[k + 1], t], t0[k] <= t < t0[k + 1]}, {k, 1, n}]];
passivationfilter = tau0/(tau0 s + 1);
Phivt[t_] := InverseLaplaceTransform[
  passivationfilter (Av s^2 + Bv s + Cv)/(s^2 + Dv s + 1) /. parms, s, t];
Phiht[t_] := InverseLaplaceTransform[
  passivationfilter (Ah s^2 + Bh s + Ch)/(s^2 + Dh s + 1) /. parms, s, t];

Xvtau = Assuming[Join[{tau, t}, U0] \[Element] Reals, 
Integrate[Phivt[tau - t]*uf[t], {t, 0, tau}]] + InverseLaplaceTransform[(s uv0 + duv0)/(s^2 + Dv s + 1) /. parms, s, tau];
Xhtau = Assuming[Join[{tau, t}, U0] \[Element] Reals, 
Integrate[Phiht[tau - t]*uf[t], {t, 0, tau}]] + InverseLaplaceTransform[(s uh0 + duh0)/(s^2 + Dh s + 1) /. parms, s, tau];

TXvtau = Table[Abs[Xvtau /. {tau -> tmin + dt/2 (k - 1)}], {k, 1, 2 n + 1}] // ComplexExpand;
TXhtau = Table[Abs[Xhtau /. {tau -> tmin + dt/2 (k - 1)}], {k, 1, 2 n + 1}] // ComplexExpand;

rho = 10;
Umax = 1;
sol = NMaximize[{TXvtau.TXvtau - rho^2 TXhtau.TXhtau, U0.U0 <= Umax}, U0, Method -> "DifferentialEvolution"]

Xvtau0 = Re[Xvtau /. sol[[2]]];
Xhtau0 = Re[Xhtau /. sol[[2]]];
Plot[{Xvtau0, Xhtau0}, {tau, tmin, tmax}, PlotStyle -> {{Blue, Thick}, {Red, Thick}}, PlotRange -> All]
U0t = Table[{tmin + dt (k - 1), u0[k]} /. sol[[2]], {k, 1, nt}];
ListLinePlot[U0t, PlotStyle -> Black, PlotRange -> All]

Follow two plots: in the first in blue $x_v$ and red $x_h$ and for the second, in black $u$.