Update optimize.py

scripts/optimize.py

@@ -1,20 +1,33 @@
 import numpy as np
 from scipy.optimize import least_squares
 
+# Эталонные значения параметров
+K1_REFERENCE = 0.577
+K2_REFERENCE = 0.422
 
-REFERENCE_K1 = 0.577
+# Диапазоны из задания
-REFERENCE_K2 = 0.422
+K1_BOUNDS = (1e-5, 1e2)
-TIMES = np.linspace(0.0, 10.0, 101)
+K2_BOUNDS = (1e-2, 1.0)
-SUBSTEPS = 10
+
+T_BEGIN = 0.0
+T_END = 10.0
+H = 0.1
+
+
+def analytical_solution(t, k1, k2):
+    """
+    Аналитическое решение задачи Коши.
+    """
+    t = np.asarray(t, dtype=float)
 
-def analytical(t, k1, k2):
     x1 = np.exp(-k1 * t)
 
-    if k1 == k2:
+    if np.isclose(k1, k2):
         x2 = k1 * t * np.exp(-k1 * t)
     else:
-        x2 = k1 / (k2 - k1) * (
+        x2 = (
-            np.exp(-k1 * t) - np.exp(-k2 * t)
+            k1 / (k2 - k1)
+            * (np.exp(-k1 * t) - np.exp(-k2 * t))
         )
 
     x3 = 1.0 - x1 - x2
@@ -22,90 +35,102 @@ def analytical(t, k1, k2):
     return np.column_stack((x1, x2, x3))
 
 
-def solve_rk4(k1, k2):
+def derivative(state, k1, k2):
+    """
+    Правая часть системы дифференциальных уравнений.
+    """
+    x1, x2, x3 = state
 
-    def rhs(x):
+    return np.array([
-        x1, x2, x3 = x
+        -k1 * x1,
+        k1 * x1 - k2 * x2,
+        k2 * x2
+    ])
 
-        return np.array([
-            -k1 * x1,
-            k1 * x1 - k2 * x2,
-            k2 * x2,
-        ])
 
-    state = np.array([1.0, 0.0, 0.0])
+def rk4_step(state, h, k1, k2):
+    """
+    Один шаг метода Рунге--Кутты четвёртого порядка.
+    """
+    q1 = derivative(state, k1, k2)
+    q2 = derivative(state + h * q1 / 2.0, k1, k2)
+    q3 = derivative(state + h * q2 / 2.0, k1, k2)
+    q4 = derivative(state + h * q3, k1, k2)
 
-    solution = np.empty((len(TIMES), 3))
+    return state + h / 6.0 * (q1 + 2.0 * q2 + 2.0 * q3 + q4)
-    solution[0] = state
 
-    for j in range(1, len(TIMES)):
-        h = (TIMES[j] - TIMES[j - 1]) / SUBSTEPS
 
-        for _ in range(SUBSTEPS):
+def rk4_solution(t, k1, k2):
-            q1 = rhs(state)
+    """
-            q2 = rhs(state + h * q1 / 2.0)
+    Решение задачи Коши методом RK4.
-            q3 = rhs(state + h * q2 / 2.0)
+    """
-            q4 = rhs(state + h * q3)
+    solution = np.zeros((len(t), 3), dtype=float)
+    solution[0] = [1.0, 0.0, 0.0]
 
-            state = state + h / 6.0 * (
+    for index in range(1, len(t)):
-                q1 + 2.0 * q2 + 2.0 * q3 + q4
+        h = t[index] - t[index - 1]
-            )
 
-        solution[j] = state
+        solution[index] = rk4_step(
+            solution[index - 1],
+            h,
+            k1,
+            k2
+        )
 
     return solution
 
 
-def residuals(parameters, reference):
+# Моменты времени из задания:
+# tj = 0.1j, j = 0, ..., 100
+time = np.linspace(T_BEGIN, T_END, 101)
+
+# Эталонные данные
+reference_solution = analytical_solution(time, K1_REFERENCE, K2_REFERENCE)
+
+
+def residuals(parameters):
+    """
+    Вектор невязок между моделью и эталонными данными.
+    """
     k1, k2 = parameters
 
-    calculated = solve_rk4(k1, k2)
+    numerical_solution = rk4_solution(time, k1, k2)
 
-    return (calculated - reference).ravel()
+    # Получается вектор из 101 * 3 = 303 невязок
+    return (numerical_solution - reference_solution).ravel()
 
 
-def main():
+# Начальное приближение
-    reference = analytical(TIMES, REFERENCE_K1, REFERENCE_K2)
+initial_parameters = np.array([0.3, 0.7])
 
-    result = least_squares(
+result = least_squares(
-        fun=residuals,
+    residuals,
-        x0=[0.3, 0.7],
+    initial_parameters,
-        args=(reference,),
+    bounds=(
-        bounds=(
+        [K1_BOUNDS[0], K2_BOUNDS[0]],
-            [1e-5, 0.01],
+        [K1_BOUNDS[1], K2_BOUNDS[1]]
-            [100.0, 1.0],
-        ),
-        method="trf",
-        max_nfev=200,
-        verbose=1,
     )
+)
 
-    k1_found, k2_found = result.x
+# result.x - solution
+found_k1, found_k2 = result.x
 
-    absolute_error_k1 = abs(k1_found - REFERENCE_K1)
+absolute_error_k1 = abs(found_k1 - K1_REFERENCE)
-    absolute_error_k2 = abs(k2_found - REFERENCE_K2)
+absolute_error_k2 = abs(found_k2 - K2_REFERENCE)
 
-    relative_error_k1 = absolute_error_k1 / REFERENCE_K1 * 100
+relative_error_k1 = absolute_error_k1 / K1_REFERENCE * 100.0
-    relative_error_k2 = absolute_error_k2 / REFERENCE_K2 * 100
+relative_error_k2 = absolute_error_k2 / K2_REFERENCE * 100.0
+criterion_value = np.sum(result.fun ** 2)
 
-    sum_squared_errors = np.sum(result.fun ** 2)
+print("Results")
-    rmse = np.sqrt(np.mean(result.fun ** 2))
+print(f"Found k1:{found_k1:.4f}")
+print(f"Theoretical k1:{K1_REFERENCE:.4f}")
+print(f"Absolute error of k1:{absolute_error_k1:.4e}")
+print(f"Relative error of k1:{relative_error_k1:.4e} %")
+print()
 
-    print(f"\nk1 = {k1_found:.4f}")
+print(f"Found k2:{found_k2:.4f}")
-    print(f"k2 = {k2_found:.4f}")
+print(f"Theoretical k2:{K2_REFERENCE:.4f}")
-
+print(f"Absolute error of k2:{absolute_error_k2:.4e}")
-    print(f"\nАбсолютное отклонение k1: {absolute_error_k1:.6e}")
+print(f"Relative error of k2:{relative_error_k2:.4e} %")
-    print(f"Абсолютное отклонение k2: {absolute_error_k2:.6e}")
+print()
-
-    print(f"\nОтносительное отклонение k1: {relative_error_k1:.6e}%")
-    print(f"Относительное отклонение k2: {relative_error_k2:.6e}%")
-
-    print(f"\nСумма квадратов отклонений: {sum_squared_errors:.6e}")
-    print(f"RMSE: {rmse:.6e}")
-
-    if not result.success:
-        print("\nПоиск завершён без подтверждённой сходимости.")
-
-
-if __name__ == "__main__":
-    main()