Проблема

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

Для этой проблемы истинные скорости реакции неизвестны. Так, необходимо наблюдать реакции и вывести уровни. Примите, что можно измерить вещества для набора времен $$t$$. От этих наблюдений соответствуйте лучшему набору скоростей реакции к измерениям.

Модель

Модель имеет шесть веществ, $$C_1$$ через $$C_6$$, это реагирует можно следующим образом:

Скорость реакции пропорциональна продукту количеств необходимых веществ. Так, если $$y_i$$ представляет количество вещества $$C_i$$, затем скорость реакции, чтобы произвести $$C_3$$ $$r_1{y}_1{y}_2$$. Точно так же скорость реакции, чтобы произвести $$C_5$$ $$r_2{y}_3{y}_4$$, и скорость реакции, чтобы произвести $$C_6$$ $$r_3{y}_3{y}_4$$.

Другими словами, дифференциальное уравнение, управляющее эволюцией системы,

Запустите дифференциальное уравнение во время 0 в точке $$y(0)=[1,1,0,1,0,0]$$. Эти начальные значения гарантируют, что все вещества реагируют полностью, вызывая $$C_1$$ через $$C_4$$ чтобы приблизиться к нулю как к, времени увеличивается.

Специальная модель в MATLAB

diffun функционируйте реализует дифференциальные уравнения в форме, готовой к решению ode45.

type diffun

function dydt = diffun(~,y,r)
dydt = zeros(6,1);
s12 = y(1)*y(2);
s34 = y(3)*y(4);

dydt(1) = -r(1)*s12;
dydt(2) = -r(1)*s12;
dydt(3) = -r(2)*s34 + r(1)*s12 - r(3)*s34;
dydt(4) = -r(2)*s34 - r(3)*s34;
dydt(5) = r(2)*s34;
dydt(6) = r(3)*s34;
end

Истинные скорости реакции $$r_1=2.5$$, $$r_2=1.2$$, и $$r_3=0.45$$. Вычислите эволюцию системы для нуля времен до пять путем вызова ode45.

rtrue = [2.5 1.2 0.45];
y0 = [1 1 0 1 0 0];
tspan = linspace(0,5);
soltrue = ode45(@(t,y)diffun(t,y,rtrue),tspan,y0);
yvalstrue = deval(soltrue,tspan);
for i = 1:6
    subplot(3,2,i)
    plot(tspan,yvalstrue(i,:))
    title(['y(',num2str(i),')'])
end

Задача оптимизации

Чтобы подготовить проблему к решению в подходе, основанном на проблеме, создайте трехэлементную переменную r оптимизации это имеет нижнюю границу 0.1 и верхняя граница 10.

r = optimvar('r',3,"LowerBound",0.1,"UpperBound",10);

Целевая функция для этой проблемы является суммой квадратов различий между решением для ОДУ параметрами r и решение истинными параметрами yvals. Чтобы выразить эту целевую функцию, сначала запишите функцию MATLAB, которая вычисляет решение для ОДУ с помощью параметров r. Этой функцией является RtoODE функция.

type RtoODE

function solpts = RtoODE(r,tspan,y0)
sol = ode45(@(t,y)diffun(t,y,r),tspan,y0);
solpts = deval(sol,tspan);
end

Использовать RtoODE в целевой функции преобразуйте функцию в выражение оптимизации при помощи fcn2optimexpr.

myfcn = fcn2optimexpr(@RtoODE,r,tspan,y0);

Выразите целевую функцию как сумму различий в квадрате между решением для ОДУ и решением истинными параметрами.

obj = sum(sum((myfcn - yvalstrue).^2));

Создайте задачу оптимизации с целевой функцией obj.

prob = optimproblem("Objective",obj);

Просмотрите проблему путем вызова show.

show(prob)

  OptimizationProblem : 

	Solve for:
       r

	minimize :
       sum(sum((RtoODE(r, extraParams{1}, extraParams{2}) - extraParams{3}).^2, 1))

       extraParams


	variable bounds:
       0.1 <= r(1) <= 10
       0.1 <= r(2) <= 10
       0.1 <= r(3) <= 10

Решите задачу

Найти параметры оптимальной подгонки r, выскажите исходное предположение r0 для решателя и вызова solve.

r0.r = [1 1 1];
[rsol,sumsq] = solve(prob,r0)

Solving problem using lsqnonlin.

Local minimum found.

Optimization completed because the size of the gradient is less than
the value of the optimality tolerance.

<stopping criteria details>

rsol = struct with fields:
    r: [3×1 double]

sumsq = 3.8669e-15

Сумма различий в квадрате является по существу нулем, означая решатель, найденный параметрами, которые вызывают решение для ОДУ совпадать с решением истинным параметрам. Так, как ожидалось решение содержит истинные параметры.

disp(rsol.r)

    2.5000
    1.2000
    0.4500

disp(rtrue)

    2.5000    1.2000    0.4500

Ограниченные наблюдения

Предположим, что вы не можете наблюдать все компоненты y, но только окончательные результаты y(5) и y(6). Можно ли получить значения всех скоростей реакции на основе этой ограниченной информации?

Чтобы узнать, измените функциональный RtoODE возвратить только пятый и шестой ODE выходные параметры. Модифицированный решатель ОДУ находится в RtoODE2.

type RtoODE2

function solpts = RtoODE2(r,tspan,y0)
solpts = RtoODE(r,tspan,y0);
solpts = solpts([5,6],:); % Just y(5) and y(6)
end

RtoODE2 функционируйте просто вызывает RtoODE и затем берет итоговые две строки выхода.

Создайте новое выражение оптимизации из RtoODE2 и переменная r оптимизации, данные об отрезке времени tspan, и начальная точка y0.

myfcn2 = fcn2optimexpr(@RtoODE2,r,tspan,y0);

Измените данные о сравнении, чтобы включать выходные параметры 5 и 6 только.

yvals2 = yvalstrue([5,6],:);

Создайте новую объективную и новую задачу оптимизации из выражения оптимизации myfcn2 и данные о сравнении yvals2.

obj2 = sum(sum((myfcn2 - yvals2).^2));
prob2 = optimproblem("Objective",obj2);

Решите задачу на основе этого ограниченного набора наблюдений.

[rsol2,sumsq2] = solve(prob2,r0)

Solving problem using lsqnonlin.

Local minimum possible.

lsqnonlin stopped because the final change in the sum of squares relative to 
its initial value is less than the value of the function tolerance.

<stopping criteria details>

rsol2 = struct with fields:
    r: [3×1 double]

sumsq2 = 2.1616e-05

Еще раз возвращенная сумма квадратов является по существу нулем. Это означает, что решатель нашел правильные скорости реакции?

disp(rsol2.r)

    1.7811
    1.5730
    0.5899

disp(rtrue)

    2.5000    1.2000    0.4500

Нет; в этом случае новые уровни очень отличаются от истинных уровней. Однако график нового решения для ОДУ по сравнению с истинными значениями показывает тот y(5) и y(6) совпадайте с истинными значениями.

figure
plot(tspan,yvals2(1,:),'b-')
hold on
ss2 = RtoODE2(rsol2.r,tspan,y0);
plot(tspan,ss2(1,:),'r--')
plot(tspan,yvals2(2,:),'c-')
plot(tspan,ss2(2,:),'m--')
legend('True y(5)','New y(5)','True y(6)','New y(6)','Location','northwest')
hold off

Чтобы идентифицировать правильные скорости реакции для этой проблемы, у вас должны быть данные для большего количества наблюдений, чем y(5) и y(6).

Постройте все компоненты решения новыми параметрами и постройте решение истинными параметрами.

figure
yvals2 = RtoODE(rsol2.r,tspan,y0);
for i = 1:6
    subplot(3,2,i)
    plot(tspan,yvalstrue(i,:),'b-',tspan,yvals2(i,:),'r--')
    legend('True','New','Location','best')
    title(['y(',num2str(i),')'])
end

Новыми параметрами, веществами $$C_1$$ и $$C_2$$ высушивайте более медленно, и вещество $$C_3$$ не накапливается так же. Но вещества $$C_4$$, $$C_5$$, и $$C_6$$ имейте точно ту же эволюцию и новыми параметрами и истинными параметрами.