ロジスティック関数とフィードバック



### 1. **目標関数（ロジスティック関数）**

$$
x_{\text{des}}(t) = \frac{K}{1 + A e^{-rt}}
$$

* 目標とする滑らかなS字曲線
* 初期値：$x_{\text{des}}(0) = \frac{K}{1 + A}$

---

### 2. **目標関数の時間微分（成長速度）**

$$
\frac{d}{dt} x_{\text{des}}(t) = r x_{\text{des}}(t) \left(1 - \frac{x_{\text{des}}(t)}{K} \right)
= \dot{x}_{\text{des}}(t)
$$

---

### 3. **非線形システムダイナミクス**

ロジスティック微分方程式にノイズと制御入力を加える：

$$
\frac{dx}{dt} = f(x) + u(t) + \eta(t)
= r x(t) \left(1 - \frac{x(t)}{K} \right) + u(t) + \eta(t)
$$

* $\eta(t)$：外乱（ノイズ）
* $u(t)$：制御入力

---

### 4. **状態フィードバック制御則**

$$
u(t) = \dot{x}_{\text{des}}(t) - f(x(t)) - K_p \cdot \left(x(t) - x_{\text{des}}(t)\right)
$$

* $\dot{x}_{\text{des}}(t)$：目標の速度
* $f(x(t))$：実際の成長速度
* $K_p$：比例フィードバックゲイン
* 誤差：$e(t) = x(t) - x_{\text{des}}(t)$

---

### 5. **数値更新式（オイラー法による離散化）**

時間刻み $\Delta t$ により：

$$
x_{k+1} = x_k + \Delta t \cdot \left( f(x_k) + u_k + \eta_k \right)
$$

---

### 6. **外乱（ノイズ）のモデル**

$$
\eta_k =
\begin{cases}
\mathcal{N}(0, \sigma^2), & 100 \leq k < 120 \\
0, & \text{otherwise}
\end{cases}
$$

---

## ✅ 全体構成まとめ

| 構成要素 | 数式 or 意味                                                                                        |
| ---- | ----------------------------------------------------------------------------------------------- |
| 目標関数 | $x_{\text{des}}(t) = \dfrac{K}{1 + A e^{-rt}}$                                                  |
| システム | $\dfrac{dx}{dt} = r x \left(1 - \dfrac{x}{K} \right) + u(t) + \eta(t)$                          |
| 制御入力 | $u(t) = \dot{x}_{\text{des}}(t) - r x \left(1 - \dfrac{x}{K} \right) - K_p(x - x_{\text{des}})$ |
| 数値更新 | $x_{k+1} = x_k + \Delta t \cdot [f(x_k) + u_k + \eta_k]$                                        |

---
```python
import numpy as np
import matplotlib.pyplot as plt

# --- パラメータ定義 ---
r = 1.0        # 成長率
K = 100        # 環境収容力
A = 9.0        # ロジスティック関数の初期条件定数
Kp = 2.0       # フィードバックゲイン
dt = 0.1       # 時間刻み
T = 20         # 総シミュレーション時間
steps = int(T / dt)
sigma = 5.0    # ノイズ標準偏差

# --- 時間軸と目標関数 ---
t = np.linspace(0, T, steps)
x_des = K / (1 + A * np.exp(-r * t))                     # 目標ロジスティック関数
dx_des = r * x_des * (1 - x_des / K)                     # 目標の時間微分

# --- 初期化 ---
x = np.zeros(steps)                                      # 実際の状態
x[0] = 10                                                # 初期値

# --- ノイズモデル η_k ---
eta = np.zeros(steps)
eta[100:120] = np.random.normal(0, sigma, size=20)       # 10s〜12sにノイズ注入

# --- メインループ：数値更新（オイラー法） ---
for k in range(steps - 1):
    f_x = r * x[k] * (1 - x[k] / K)                      # f(x_k)
    error = x[k] - x_des[k]                              # e_k = x_k - x_des_k
    u = dx_des[k] - f_x - Kp * error                     # u_k = dx_des - f(x) - Kp * e
    x[k+1] = x[k] + dt * (f_x + u + eta[k])              # x_{k+1} = x_k + dt * [f + u + η]

# --- プロット ---
plt.figure(figsize=(10, 5))
plt.plot(t, x, label="Controlled $x(t)$ with Noise", linewidth=2)
plt.plot(t, x_des, '--', label="Target $x_{\\mathrm{des}}(t)$", linewidth=2)
plt.axvspan(10, 12, color='gray', alpha=0.2, label="Noise Zone")
plt.xlabel("Time")
plt.ylabel("Population $x(t)$")
plt.title("Tracking Logistic Function with Feedback and Noise Rejection")
plt.grid(True)
plt.legend()
plt.tight_layout()
plt.show()

```
結果
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/4046871/d70bb53d-7c42-4772-ac89-e3e47b9491fc.png)