I am looking for references on the implementation of computational optimization methods within the context of differential equations.
As an example, let's say I have several coupled non-linear ordinary differential equations describing a biogeochemical system where some living organism $X$ eats chemical compounds $c_1$, $c_2$ and $c_3$ that contain some energy for living and reproduction ($\Omega_1$, $\Omega_2$ and $\Omega_3$). This energy is variable in time and depends on the concentrations of $c_1$, $c_2$ and $c_3$:
$$\begin{aligned} \dot c_1 &= - X \cdot c_1 \cdot \alpha\\ \dot c_2 &= - X \cdot c_2 \cdot \beta \\ \dot c_3 &= - X \cdot c_3 \cdot \gamma\\ \dot X &= \Omega_1 (c_1, c_2 ,c_3)\cdot X \cdot c_1 \cdot \alpha + \Omega_2(c_1, c_2 ,c_3)\cdot X \cdot c_2 \cdot \beta + \Omega_3(c_1, c_2 ,c_3)\cdot X \cdot c_3 \cdot \gamma \end{aligned}$$
where
$$\alpha + \beta + \gamma = 1$$
At every timestep, an optimization algorithm finds the maximum
$$\max \bigg( \Omega_1 (c_1)\cdot X \cdot c_1 \cdot \alpha + \Omega_2(c_2)\cdot X \cdot c_2 \cdot \beta + \Omega_3(c_3)\cdot X \cdot c_3 \cdot \gamma \bigg)$$
i.e., to find the best value of $\alpha, \beta, \gamma$ to ensure maximum growth of $X$.
The problem I am facing is that the solution is not smooth, i.e., on some time step living organism picks best energy-yielding product $c_1$, but the next time step it picks $c_2$, so it produces non-smooth solutions for $\alpha, \beta, \gamma$ and it kind of jumps back and forth between $c_1$, $c_2$ and $c_3$ in time, which is not very realistic. So I need some kind of optimization algorithms that has smooth switching functions, e.g., sigmoid functions, but for multiple variables $\alpha, \beta$, and $\gamma$.
This is a toy example. In the real example, I am solving PDE of nonlinear equations and have more than $30$ variables of $\alpha$, $\beta$, and $\gamma$.
I found some information about receding horizon algorithm for the system that is very close to the system I described above (link). However, the description is very general and no information on the implementation.
Hint: I will assume that $\alpha$, $\beta$ and $\gamma$ are constants. You can simplify your problem by eliminating $c_2$ and $c_3$ by dividing the first and second and the first and the third differential equation.
$$\dfrac{dc_1}{dc_2}=\dfrac{\alpha}{\beta}\dfrac{c_1}{c_2}\implies \dfrac{dc_1}{c_1}=\dfrac{\alpha}{\beta}\dfrac{dc_2}{c_2}$$ $$\implies \ln c_1 -\ln c_1(t=0)=\dfrac{\alpha}{\beta}[\ln c_2 -\ln c_2(t=0)]$$ Now, solve for $c_2$ to obtain
$$c_2=c_2(0)\left[\frac{c_1}{c_1(0)} \right]^{\beta/\alpha}.$$
Do the same procedure for the first and third equation. If the parameters are not constant then you will have to solve a time-variant system. Then you will only be left with a second order system. You should be able to use a similar method for the resulting system.