Проблемы с интерполяцией и поиском минимального значения

[Mr.Audar]

Всем доброго времени суток. Задача стоит в построении математической модели полимеризации полиизопрена и нахождения оптимального начального значения ДИБАГ (т.е s[:, 2]), который даст Среднемассовое значение (Mw) = 600000.
В теории начальное значение s[:, 2] (концентрация реагента) должна быть около 0,00378, чтобы получить Mw = 60000. На данном этапе, мой код, если его можно так назвать не показывает мне должного результата. По поиску минимального значения мне подсказали, вставил как в инструкции. Прошу помощи данного форума. Представляю сам код и вывод.

Используется синтаксис Python

import matplotlib.pyplot as plt
import numpy as np
from scipy.optimize import minimize
from scipy.integrate import odeint
from scipy.interpolate import splev, splrep, sproot
 
def  model_polisopren (s,t): # в данной функции записана Сист.Диф.Ур хим.процесса
    kp = 48; km = 0.0048; ka1 = 8.16; ka2 = 0.96 #константы скорости реакции
    dMdt = -s[0]*s[3]*(kp+km)-s[0]*(kp+km)*s[5]
    dA1dt = -ka1*s[1]*s[3]-ka1*s[1]*s[5]
    dA2dt = -ka2*s[2]*s[3]-ka2*s[2]*s[5]
    dPdt = -kp*s[0]*s[3]+(km*s[0]+ka1*s[1]+ka2*s[2])*s[5]
    dQdt = km*s[0]*s[3]+(ka1*s[1]+ka2*s[2])*s[3]
    dm0dt = kp*s[0]*s[3]-(km*s[0]+(ka1*s[1]+ka2*s[2]))*s[5]
    dn0dt = (km*s[0]+ka1*s[1]+ka2*s[2])*s[5]
    dm1dt = 2*kp*s[0]*s[3]+kp*s[0]*s[5]-(km*s[0]+(ka1*s[1]+ka2*s[2]))*s[7]
    dn1dt = (km*s[0]+(ka1*s[1]+ka2*s[2]))*s[7]
    dm2dt = 4*kp*s[0]*s[3]+kp*s[0]*s[5]+2*kp*s[0]*s[7]-(km*s[0]+(ka1*s[1]+ka2*s[2]))*s[9]
    dn2dt = (km*s[0]+(ka1*s[1]+ka2*s[2]))*s[9]
    return [dMdt, dA1dt, dA2dt, dPdt, dQdt, dm0dt, dn0dt, dm1dt, dn1dt, dm2dt, dn2dt]
    
s0 = [1.39,0.000177,0.00168,0.000007,0,0,0,0,0,0,0] # начальные значение концентрации реагентов
t = np.linspace(0, 5100, 10000)                     # задается время и шаг
s = odeint(model_polisopren, s0, t)                 # численное решение СДУ
Mw = ((s[:,9]+s[:,10])/(s[:,7]+s[:,8]+2**-100))*68  # считаем Mw тут 
plt.plot(t,Mw,'b-',linewidth=2.0,label='Mw')
plt.xlabel("t")
plt.ylabel("Mw")
plt.legend()
plt.grid()
plt.show()
    
def poisk_min (s0):  # данная функция должна находить значения нач.концетрации при определенных значений
    x = splrep(t, Mw - 600000) # задаем условия, где нач.конц. должны быть такими, чтобы Mw = 60000 (мы считали его ранее)
    x0=sproot(x)
    y1=splrep(t, s[:, 2]) # указываем какая нач.концентрация нас интересует 
    s1=splev(x0[0], y1)
    print (s1)
    return s1
res = minimize(poisk_min, s0, method='Nelder-Mead', tol=1e-3) 
print (res)

И выводит следующие значения:

0.001676790666781049
 final_simplex: (array([[1.39000000e+00, 1.77000000e-04, 1.68000000e-03, 7.00000000e-06,
        0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,
        0.00000000e+00, 0.00000000e+00, 0.00000000e+00],
       [1.39054297e+00, 1.77000000e-04, 1.68000000e-03, 7.00000000e-06,
        0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,
        0.00000000e+00, 0.00000000e+00, 0.00000000e+00],
       [1.39000000e+00, 1.77069141e-04, 1.68000000e-03, 7.00000000e-06,
        0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,
        0.00000000e+00, 0.00000000e+00, 0.00000000e+00],
       [1.39000000e+00, 1.77000000e-04, 1.68065625e-03, 7.00000000e-06,
        0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,
        0.00000000e+00, 0.00000000e+00, 0.00000000e+00],
       [1.39000000e+00, 1.77000000e-04, 1.68000000e-03, 7.00273437e-06,
        0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,
        0.00000000e+00, 0.00000000e+00, 0.00000000e+00],
       [1.39000000e+00, 1.77000000e-04, 1.68000000e-03, 7.00000000e-06,
        1.95312500e-06, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,
        0.00000000e+00, 0.00000000e+00, 0.00000000e+00],
       [1.39000000e+00, 1.77000000e-04, 1.68000000e-03, 7.00000000e-06,
        0.00000000e+00, 1.95312500e-06, 0.00000000e+00, 0.00000000e+00,
        0.00000000e+00, 0.00000000e+00, 0.00000000e+00],
       [1.39000000e+00, 1.77000000e-04, 1.68000000e-03, 7.00000000e-06,
        0.00000000e+00, 0.00000000e+00, 1.95312500e-06, 0.00000000e+00,
        0.00000000e+00, 0.00000000e+00, 0.00000000e+00],
       [1.39000000e+00, 1.77000000e-04, 1.68000000e-03, 7.00000000e-06,
        0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 1.95312500e-06,
        0.00000000e+00, 0.00000000e+00, 0.00000000e+00],
       [1.39000000e+00, 1.77000000e-04, 1.68000000e-03, 7.00000000e-06,
        0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,
        1.95312500e-06, 0.00000000e+00, 0.00000000e+00],
       [1.39000000e+00, 1.77000000e-04, 1.68000000e-03, 7.00000000e-06,
        0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,
        0.00000000e+00, 1.95312500e-06, 0.00000000e+00],
       [1.39000000e+00, 1.77000000e-04, 1.68000000e-03, 7.00000000e-06,
        0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,
        0.00000000e+00, 0.00000000e+00, 1.95312500e-06]]), array([0.00167679, 0.00167679, 0.00167679, 0.00167679, 0.00167679,
       0.00167679, 0.00167679, 0.00167679, 0.00167679, 0.00167679,
       0.00167679, 0.00167679]))
           fun: 0.001676790666781049
       message: 'Optimization terminated successfully.'
          nfev: 103
           nit: 8
        status: 0
       success: True
             x: array([1.39e+00, 1.77e-04, 1.68e-03, 7.00e-06, 0.00e+00, 0.00e+00,
       0.00e+00, 0.00e+00, 0.00e+00, 0.00e+00, 0.00e+00])
​

[kotenok gav]

А вы уверены, что то, что ваша функция model_polisopren не использует t - правильно?

[Mr.Audar]

kotenok gav, да аргумент t используется правильно. В нем задается ось X с размерностью от 0 до 5100 с шагом 10000. (5100 это время, т.е секунды)

[kotenok gav]

Я не об этом. Он нигде не используется в теле функции.

[Mr.Audar]

kotenok gav, Аргумент t задается в Odeint для численного решения и затем для построения графика ))))

[12d3]

Кажется, ваша минимизация делает совсем не то, что нужно. Уточните, требуется подобрать такое начальное значение s0[2] при неизменных остальных начальных значениях, чтобы Mw через 5100 секунд было как можно ближе к 600000?

[Mr.Audar]

12d3
Да, моя минимизация немного неправильная ((( А так, задача верная - нужно найти s0[2] при неизменных начальных значениях.

[12d3]

Тогда, по сути, вам нужно численно решить уравнение $f(x) = 0$, где $f$ такова:

Используется синтаксис Python

def f(x):
    s0 = [1.39, 0.000177, 0.00168, 0.000007, 0, 0, 0, 0, 0, 0, 0]
    s0[2] = x
    t = np.linspace(0, 5100, 10000)               
    s = odeint(model_polisopren, s0, t)                 
    Mw = ((s[:,9]+s[:,10])/(s[:,7]+s[:,8]+2**-100))*68
    return Mw[-1] - 600000
 

Для этого есть подходящая функция из пакета scypi: https://docs.scipy.org/doc/scipy/reference/generated/scipy.optimize.fsolve.html

[Mr.Audar]

12d3
Спасибо. Буду изучать )))) Вы не подскажите какой метод можно реализовать для двумерной минимизации ??