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

@Koala1996 · Регистрация: 04.12.2019

Студворк — интернет-сервис помощи студентам

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

Python
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
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)

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

Python
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
0.0037687965538383606
 final_simplex: (array([[1.39000000e+00, 1.77000000e-04, 3.78000000e-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, 3.78000000e-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, 3.78000000e-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, 3.78147656e-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, 3.78000000e-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, 3.78000000e-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, 3.78000000e-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, 3.78000000e-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, 3.78000000e-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, 3.78000000e-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, 3.78000000e-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, 3.78000000e-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.0037688, 0.0037688, 0.0037688, 0.0037688, 0.0037688, 0.0037688,
       0.0037688, 0.0037688, 0.0037688, 0.0037688, 0.0037688, 0.0037688]))
           fun: 0.0037687965538383606
       message: 'Optimization terminated successfully.'
          nfev: 103
           nit: 8
        status: 0
       success: True
             x: array([1.39e+00, 1.77e-04, 3.78e-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])
​

@u235 · 24.11.2021, 19:20

Химическая кинетика, это интересно, в свое время проходил курс.

А у вас точно с размерностями констант скорости и концентраций все в порядке? В чем они у вас измеряются? Просто их различие на порядки как-то странно.

@Koala1996 · 24.11.2021, 19:39 **[ТС]**

Да, константы скоростей верные. За основу была взята модель со всеми вытекающими значениями из научной статьи (прилагаю во вложение). Тут реализовали с помощью генетического метода. Мне поставили задачу использовать математический метод ))))

@Koala1996 · 24.11.2021, 20:30 **[ТС]**

u235, приятно осознавать, что нашел коллегу на этом форуме

У вас есть какие-нибудь идеи по поводу решения ??? Мне также предлагали реализовать метод стрельбы.

@u235 · 24.11.2021, 21:08

Думаю что должно в итоге сойтись с методом, описанным в статье..
Koala1996, не совсем понятно, в статье 13 уравнений, а у вас в скрипте 11... где еще 2 уравнения?
С начальными условиями более-менее понятно, а что с константами скоростей? Откуда из взяли?

@Koala1996 · 24.11.2021, 21:24 **[ТС]**

u235, да, все верно, генетический метод может выполнить данную задачу. И на сколько мне известно он не стабилен и узкоспециализирован под определенные хим.реакции. Моя задача сделать некоторый универсальный скрипт полимеризационных процессов. Почему не 13, а 11, потому что оставшиеся 2 для других задач.
Авторы данной статьи сказали о значений константы и сказали принять за истину ))))

@u235 · 24.11.2021, 23:26

Сообщение от Koala1996

(концентрация реагента) должна быть около 0,00378, чтобы получить Mw = 60000.

не 0.00378, а 0.000378, и да, тогда ваш скрипт (график если быть точнее) выдает молекулярный вес 60000 за 85 минут, так что константы правильные и я зря сомневался.
Теперь надо разобраться с поиском минимума. Получается что нужно минимизировать функцию одной переменной - начальная концентрация ДИБАГ s2[0]? Зачем использовать симплекс-метод, если переменная одна..?

Добавлено через 25 минут
Вот, выдает что-то близкое к 0.000378:

Python
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
import matplotlib.pyplot as plt
import numpy as np
from scipy.optimize import  minimize_scalar
from scipy.integrate import odeint
 
 
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]
    
def funct1(s2):
    s0 = [1.39,s2,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[-1,9]+s[-1,10])/(s[-1,7]+s[-1,8]+2**-100))*68  # считаем Mw тут 
    return np.abs(Mw-Mw_set)
    
Mw_set=600000
 
res = minimize_scalar(funct1, bounds=(0.0001, 0.001), method='bounded', tol=1e-6) 
print (res)

@Koala1996 · 25.11.2021, 09:23 **[ТС]**

Симплекс метод использовал аргументировав это тем, что используется итерационной метод. Перечитав теоретическую базу, я понял, что скорее всего был не прав.
u235, большое спасибо за помощь ! Буду изучать ваш код, т.к некоторые моменты для меня еще непонятны.

Новые блоги и статьи Все статьи Все блоги /
Программный контроль заполнения реквизита табличной части документа Maks 02.04.2026 Алгоритм из решения ниже реализован на примере нетипового документа "СписаниеМатериалов", разработанного в конфигурации КА2. Задача: реализовать контроль заполнения реквизита "ПричинаСписания". . .	wmic не является внутренней или внешней командой Maks 02.04.2026 Решение: DISM / Online / Add-Capability / CapabilityName:WMIC~~~~ Отсюда: https:/ / winitpro. ru/ index. php/ 2025/ 02/ 14/ komanda-wmic-ne-naydena/	Программная установка даты и запрет ее изменения Maks 02.04.2026 Алгоритм из решения ниже реализован на примере нетипового документа "СписаниеМатериалов", разработанного в конфигурации КА2. Задача: при создании документов установить период списания автоматически. . .	Вывод данных в справочнике через динамический список Maks 01.04.2026 Реализация из решения ниже выполнена на примере нетипового справочника "Спецтехника" разработанного в конфигурации КА2. Задача: вывести данные из ТЧ нетипового документа. . .
Программное заполнения текстового поля в реквизите формы документа Maks 01.04.2026 Алгоритм из решения ниже реализован на нетиповом документе "ВыдачаОборудованияНаСпецтехнику" разработанного в конфигурации КА2, в дополнении к предыдущему решению. На форме документа создается. . .	К слову об оптимизации kumehtar 01.04.2026 Вспоминаю начало 2000-х, университет, когда я писал на Delphi. Тогда среди программистов на форумах активно обсуждали аккуратную работу с памятью: нужно было следить за переменными, вовремя. . .	Идея фильтра интернета (сервер = слой+фильтр). Hrethgir 31.03.2026 Суть идеи заключается в том, чтобы запустить свой сервер, о чём я если честно мечтал давно и давно приобрёл книгу как это сделать. Но не было причин его запускать. Очумелые учёные напечатали на. . .	Модель здравосоХранения 6. ESG-повестка и устойчивое развитие; углублённый анализ кадрового бренда anaschu 31.03.2026 В прикрепленном документе раздумья о том, как можно поменять модель в будущем

@u235 5516 / 2869 / 571 Регистрация: 07.11.2019 Сообщений: 4,760
	24.11.2021, 19:20
	Химическая кинетика, это интересно, в свое время проходил курс. А у вас точно с размерностями констант скорости и концентраций все в порядке? В чем они у вас измеряются? Просто их различие на порядки как-то странно. 0

@Koala1996 0 / 0 / 0 Регистрация: 04.12.2019 Сообщений: 15
	24.11.2021, 20:30 [ТС]
	u235, приятно осознавать, что нашел коллегу на этом форуме У вас есть какие-нибудь идеи по поводу решения ??? Мне также предлагали реализовать метод стрельбы. 0

@u235 5516 / 2869 / 571 Регистрация: 07.11.2019 Сообщений: 4,760
	24.11.2021, 21:08
	Думаю что должно в итоге сойтись с методом, описанным в статье.. Koala1996, не совсем понятно, в статье 13 уравнений, а у вас в скрипте 11... где еще 2 уравнения? С начальными условиями более-менее понятно, а что с константами скоростей? Откуда из взяли? 0

@Koala1996 0 / 0 / 0 Регистрация: 04.12.2019 Сообщений: 15
	24.11.2021, 21:24 [ТС]
	u235, да, все верно, генетический метод может выполнить данную задачу. И на сколько мне известно он не стабилен и узкоспециализирован под определенные хим.реакции. Моя задача сделать некоторый универсальный скрипт полимеризационных процессов. Почему не 13, а 11, потому что оставшиеся 2 для других задач. Авторы данной статьи сказали о значений константы и сказали принять за истину )))) 0

@Koala1996 0 / 0 / 0 Регистрация: 04.12.2019 Сообщений: 15
	25.11.2021, 09:23 [ТС]
	Симплекс метод использовал аргументировав это тем, что используется итерационной метод. Перечитав теоретическую базу, я понял, что скорее всего был не прав. u235, большое спасибо за помощь ! Буду изучать ваш код, т.к некоторые моменты для меня еще непонятны. 0

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

Решение