seimeikagaku2-1.py

# https://opensource.org/licenses/GPL-3.0
import control.optimal as opt
import control as ct
import numpy as np
import pandas as pd
from numpy.typing import NDArray

X0 = np.array([1, 1])
Tf = 1
horizon = np.linspace(0, Tf, 15, endpoint=True)
params = {'r': 3}


def model_update(t:float, x:"NDArray", u:"NDArray", params: dict) -> "NDArray":
    """Returns dx/dt = f(x, u, t).
    """
    r = params['r']
    u[0] = np.clip(u[0], 0, 1)
    return np.array([
        (1-u[0])*x[0]*r,
        u[0]*x[0]*r,
    ])


def model_output(t:float, x:"NDArray", u:"NDArray", params: dict) -> "NDArray":
    """Returns x."""
    return x


model = ct.NonlinearIOSystem(model_update, model_output, params=params,
                             inputs=('u',), outputs=('C', 'R'), states=2)
constraints = [opt.input_range_constraint(model, [0., ], [1., ])]

def cost(x:"NDArray", u:"NDArray") -> float:
    """Stage cost function."""
    return 0

def term_cost(x:"NDArray", u:"NDArray") -> float:
    """Terminal cost function."""
    return -x[1]

result = opt.solve_ocp(model, terminal_cost = term_cost, cost=cost,
                       constraints=constraints, X0=X0, horizon=horizon)

u = result.inputs
t, X = ct.input_output_response(model, horizon, u, X0)

# plot
df = pd.DataFrame(np.concatenate([u, X]).T, index=t, columns=['u', 'C', 'R'])
df.plot()