Химический реактор с питанием

В периодическом реакторе [1] происходят следующие необратимые и экзотермические реакции:

A + B => C  (desired product)
C => D      (undesired product)

Партия начинается с реактора, частично заполненного известными концентрациями реагентов A и B. Партия реагирует в течение 0,5 часов, в течение которых дополнительные B может быть добавлена и температура реактора может быть изменена.

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

Манипулированные переменные

Измеренное нарушение порядка

Государства

Выходы

Цель состоит в том, чтобы максимизировать производство C (y1) в конце пакетного процесса. В процессе пакетной обработки должны быть удовлетворены следующие рабочие ограничения:

Задайте номинальное рабочее условие в начале пакетного процесса.

c_A0 = 10;
c_B0 = 1.167;
c_C0 = 0;
V0  = 1;
T0 = 50 + 273.15;
c_Bin = 20;

Задайте номинальные состояния.

x0 = zeros(3,1);
x0(1) = c_A0*V0;
x0(2) = x0(1) + c_C0*V0;
x0(3) = V0;
x0(4) = T0;

Задайте номинальные входы.

u0 = zeros(3,1);
u0(2) = 40;
u0(3) = c_Bin;

Задайте номинальные выходы.

y0 = fedbatch_OutputFcn(x0,u0);

Нелинейный проект MPC для оптимизации пакетной операции

Создайте нелинейный объект MPC с 4 состояниями, 3 выходами, 2 манипулированными переменными и 1 измеренным нарушением порядка.

nlmpcobj_Plan = nlmpc(4, 3, 'MV', [1,2], 'MD', 3);

In standard cost function, zero weights are applied by default to one or more OVs because there are fewer MVs than OVs.

Учитывая ожидаемую длительность пакета Tf, выберите значение шага расчета контроллера Ts и горизонт предсказания.

Tf = 0.5;
N = 50;
Ts = Tf/N;
nlmpcobj_Plan.Ts = Ts;
nlmpcobj_Plan.PredictionHorizon = N;

Если вы задаете контрольный горизонт равным горизонту предсказания, будет 50 свободных движений управления, что приводит в общей сложности к 100 переменным решения, потому что у объекта есть две манипулированные переменные. Чтобы уменьшить количество переменных принятия решений, можно задать горизонт управления с помощью блокирующих перемещений. Разделите горизонт предсказания на 8 блоков, что представляет собой 8 свободных ходов управления. Каждый из первых семи блоков длится семь шагов предсказания. Это уменьшает количество переменных принятия решений до 16.

nlmpcobj_Plan.ControlHorizon = [7 7 7 7 7 7 7 1];

Задайте нелинейную модель в контроллере. Функция fedbatch_StateFcnDT преобразует модель непрерывного времени в дискретное время с помощью многоэтапной формулы интегрирования Forward Euler.

nlmpcobj_Plan.Model.StateFcn = @(x,u) fedbatch_StateFcnDT(x,u,Ts);
nlmpcobj_Plan.Model.OutputFcn = @(x,u) fedbatch_OutputFcn(x,u);
nlmpcobj_Plan.Model.IsContinuousTime = false;

Задайте ограничения для скорости подачи B.

nlmpcobj_Plan.MV(1).Min = 0;
nlmpcobj_Plan.MV(1).Max = 1;

Задайте границы уставки температуры реактора.

nlmpcobj_Plan.MV(2).Min = 20;
nlmpcobj_Plan.MV(2).Max = 50;

Задайте верхнюю границу для скорости отвода тепла. Истинное ограничение 1.5e5. Поскольку нелинейный MPC может применять ограничения только в моменты выборки, используйте запас прочности 0.05e5 чтобы предотвратить нарушение ограничений между моментами дискретизации.

nlmpcobj_Plan.OV(2).Max = 1.45e5;

Задайте верхнюю границу для объема жидкости в реакторе.

nlmpcobj_Plan.OV(3).Max = 1.1;

Поскольку цель состоит в том, чтобы максимизировать y1, количество C в реакторе в конце пакетного времени задайте пользовательскую функцию затрат, которая заменяет квадратичные затраты по умолчанию. Начиная с y1 = x2-x1, определите пользовательскую стоимость, которая будет минимизирована следующим образом x1-x2 использование анонимной функции.

nlmpcobj_Plan.Optimization.CustomCostFcn = @(X,U,e,data) X(end,1)-X(end,2);
nlmpcobj_Plan.Optimization.ReplaceStandardCost = true;

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

nlmpcobj_Plan.Optimization.MVInterpolationOrder = 1;

Используйте нелинейное программирование по умолчанию решателя fmincon для решения нелинейной задачи MPC. В этом примере установите допуск шага решателя, чтобы помочь достичь оптимальности первого порядка.

nlmpcobj_Plan.Optimization.SolverOptions.StepTolerance = 1e-8;

Перед проведением оптимизации проверьте, удовлетворяют ли все пользовательские функции требованиям NLMPC, используя validateFcns команда.

validateFcns(nlmpcobj_Plan, x0, u0(1:2), u0(3));

Model.StateFcn is OK.
Model.OutputFcn is OK.
Optimization.CustomCostFcn is OK.
Analysis of user-provided model, cost, and constraint functions complete.

Анализ результатов оптимизации

Найдите оптимальные траектории для манипулируемых переменных, такие что производство C максимизируется в конце пакетного процесса. Для этого используйте nlmpcmove функция.

fprintf('\nOptimization started...\n');
[~,~,Info] = nlmpcmove(nlmpcobj_Plan,x0,u0(1:2),zeros(1,3),u0(3));
fprintf('   Expected production of C (y1) is %g moles.\n',Info.Yopt(end,1));
fprintf('   First order optimality is satisfied (Info.ExitFlag = %i).\n',...
    Info.ExitFlag);
fprintf('Optimization finished...\n');

Optimization started...
Slack variable unused or zero-weighted in your custom cost function. All constraints will be hard.
   Expected production of C (y1) is 2.02353 moles.
   First order optimality is satisfied (Info.ExitFlag = 1).
Optimization finished...

Дискретизированная модель использует простое интегрирование Эйлера, которое может быть неточным. Чтобы проверить это, интегрируйте модель с помощью ode15s команда для вычисленной оптимальной траектории СН.

Nstep = size(Info.Xopt,1) - 1;
t = 0;
X = x0';
t0 = 0;
for i = 1:Nstep
    u_in = [Info.MVopt(i,1:2)'; c_Bin];
    ODEFUN = @(t,x) fedbatch_StateFcn(x, u_in);
    TSPAN = [t0, t0+Ts];
    Y0 = X(end,:)';
    [TOUT,YOUT] = ode15s(ODEFUN,TSPAN,Y0);
    t = [t; TOUT(2:end)];
    X = [X; YOUT(2:end,:)];
    t0 = t0 + Ts;
end
nx = size(X,1);
Y = zeros(nx,3);
for i = 1:nx
    Y(i,:) = fedbatch_OutputFcn(X(i,:)',u_in)';
end
fprintf('\n   Actual Production of C (y1) is %g moles.\n',X(end,2)-X(end,1));
fprintf('   Heat removal rate (y2) satisfies the upper bound.\n');

   Actual Production of C (y1) is 2.0228 moles.
   Heat removal rate (y2) satisfies the upper bound.

На верхнем графике следующей рисунка, фактическое производство C согласен с ожидаемым производством C вычисляется из nlmpcmove. На нижнем графике скорость отвода тепла никогда не превышает своего жесткого ограничения.

figure
subplot(2,1,1)
plot(t,Y(:,1),(0:Nstep)*Ts, Info.Yopt(:,1),'*')
axis([0 0.5 0 Y(end,1) + 0.1])
legend({'Actual','Expected'},'location','northwest')
title('Mol C in reactor (y1)')
subplot(2,1,2)
tTs = (0:Nstep)*Ts;
t(end) = 0.5;
plot(t,Y(:,2),'-',[0 tTs(end)],1.5e5*ones(1,2),'r--')
axis([0 0.5 0.8e5, 1.6e5])
legend({'q_r','Upper Bound'},'location','southwest')
title('Heat removal rate (y2)')

Тщательный анализ скорости отвода тепла показывает, что он может проявлять peaks и овраги между моментами отбора проб, когда изменяются композиции реагентов. Следовательно, скорость отвода тепла превышает заданный максимум 1.45e5 (около t = 0,35 ч), но остается ниже истинного максимума 1.5e5.

Следующий рисунок показывает оптимальную траекторию запланированных корректировок в B скорость подачи (u1), и температура реактора (x4) и его уставку (u2).

figure
subplot(2,1,1)
stairs(tTs,Info.MVopt(:,1))
title('Feed rate of B (u1)')
subplot(2,1,2)
plot(tTs,Info.MVopt(:,2),'*',t,X(:,4)-273.15,'-',...
    [0 0.5],[20 20],'r--',[0 0.5],[50 50],'r--')
axis([0 0.5 15 55])
title('Reactor temperature and its setpoint')
legend({'Setpoint','Actual'},'location','southeast')

Траектория начинается с относительно высокой скорости подачи, которая увеличивается c_B и результирующие C скорость производства. Чтобы предотвратить превышение скорости отвода тепла, температура реактора и скорость подачи должны уменьшиться. Температура в конечном счете достигает своей нижней границы и остается там, пока реактор не будет почти полным и B скорость подачи должна быть равной нулю. Температура затем увеличивается до максимума (чтобы увеличить C производство) и, наконец, немного падает (для уменьшения производства D, что благоприятствуется при более высоких температурах).

Верхний график следующего рисунка показывает потребление c_A, что имеет тенденцию к снижению C производство. Чтобы компенсировать, план сначала увеличивается c_Bи когда это больше невозможно (объем жидкости реактора не должен превышать 1,1), в плане оптимально используется температура. На нижнем графике следующего рисунка объем жидкости никогда не превышает верхнюю границу.

figure
subplot(2,1,1)
c_A = X(:,1)./X(:,3);
c_B = (c_Bin*X(:,3) + X(:,1) + V0*(c_B0 - c_A0 - c_Bin))./X(:,3);
plot(t,[c_A, c_B])
legend({'c_A','c_B'}, 'location', 'west')
subplot(2,1,2)
plot(tTs,Info.Yopt(:,3))
title('Liquid volume')

Нелинейный проект MPC для отслеживания оптимальной траектории продукта C

Чтобы отследить оптимальную траекторию продукта C Вычисленный выше, вы проектируете другой нелинейный контроллер MPC с той же моделью предсказания и ограничениями. Однако используйте стандартные квадратичные затраты и горизонты по умолчанию в целях отслеживания.

Чтобы упростить задачу управления, примите, что оптимальная траектория B Скорость подачи реализована в объекте, и контроллер отслеживания рассматривает ее как измеренное нарушение порядка. Поэтому контроллер использует уставку температуры реактора в качестве своей единственной манипулированной переменной, чтобы отслеживать желаемую y1 профиль.

Создайте контроллер отслеживания.

nlmpcobj_Tracking = nlmpc(4,3,'MV',2,'MD',[1,3]);
nlmpcobj_Tracking.Ts = Ts;
nlmpcobj_Tracking.Model = nlmpcobj_Plan.Model;
nlmpcobj_Tracking.MV = nlmpcobj_Plan.MV(2);
nlmpcobj_Tracking.OV = nlmpcobj_Plan.OV;
nlmpcobj_Tracking.Weights.OutputVariables = [1 0 0];        % track y1 only
nlmpcobj_Tracking.Weights.ManipulatedVariablesRate = 1e-6;  % agressive MV

In standard cost function, zero weights are applied by default to one or more OVs because there are fewer MVs than OVs.

Получите C производство (y1) опорный сигнал от оптимальной траектории плана.

Cref = Info.Yopt(:,1);

Получите скорость подачи B (u1) от оптимальной траектории плана. Концентрация сырья B (u3) является константой.

MD = [Info.MVopt(:,1) c_Bin*ones(N+1,1)];

Сначала запустите контроллер слежения в нелинейном режиме MPC.

[X1,Y1,MV1,et1] = fedbatch_Track(nlmpcobj_Tracking,x0,u0(2),N,Cref,MD);
fprintf('\nNonlinear MPC: Elapsed time = %g sec. Production of C = %g mol\n',et1,Y1(end,1));

Nonlinear MPC: Elapsed time = 2.94296 sec. Production of C = 2.01131 mol

Во-вторых, запустите контроллер как адаптивный контроллер MPC.

nlmpcobj_Tracking.Optimization.RunAsLinearMPC = 'Adaptive';
[X2,Y2,MV2,et2] = fedbatch_Track(nlmpcobj_Tracking,x0,u0(2),N,Cref,MD);
fprintf('\nAdaptive MPC: Elapsed time = %g sec. Production of C = %g mol\n',et2,Y2(end,1));

Adaptive MPC: Elapsed time = 1.67153 sec. Production of C = 2.01567 mol

В-третьих, запустите контроллер как изменяющийся во времени контроллер MPC.

nlmpcobj_Tracking.Optimization.RunAsLinearMPC = 'TimeVarying';
[X3,Y3,MV3,et3] = fedbatch_Track(nlmpcobj_Tracking,x0,u0(2),N,Cref,MD);
fprintf('\nTime-varying MPC: Elapsed time = %g sec. Production of C = %g mol\n',et3,Y3(end,1));

Time-varying MPC: Elapsed time = 1.46781 sec. Production of C = 2.02349 mol

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

В этом примере все три метода приближаются к оптимальному C производство, полученное на стадии планирования.

figure
plot(Ts*(0:N),[Y1(:,1) Y2(:,1) Y3(:,1)])
title('Production of C')
legend({'NLMPC','Adaptive','TimeVarying'},'location','northwest')

Неожиданным результатом является то, что изменяющийся во времени MPC создает больше C чем нелинейный MPC. Объяснение состоит в том, что методы линеаризации модели, используемые в адаптивном и изменяющемся во времени режимах, приводят к нарушению ограничения отвода тепла, что приводит к более высокой C производство.

figure
plot(Ts*(0:N),[Y1(:,2) Y2(:,2) Y3(:,2) 1.5e5*ones(N+1,1)])
title('Heat removal rate')
legend({'NLMPC','Adaptive','TimeVarying','Constraint'},'location','southwest')

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

Изменяющийся во времени метод избегает этой проблемы. Однако в начале партии он принимает (по умолчанию), что состояния останутся постоянными на горизонте. Он исправляет для этого, когда получает свое первое решение (используя данные в opts переменная), но его первоначальный выбор температуры реактора слишком высок, что приводит к раннему q_r нарушение ограничений.

Ссылки

[1] Srinivasan, B., S. Palanki, and D. Bonvin, «Динамическая оптимизация пакетных процессов I. Характеристика номинального решения», Computers and Chemical Engineering, vol. 27 (2003), pp. 1-26.