Питаемый химический реактор периодического действия

В реакторе периодического действия [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 преобразует модель непрерывного времени в дискретное время, используя многошаговую формулу интеграции Эйлера.

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)')

Тщательное изучение скорости отвода тепла показывает, что он может иметь пики и впадины между моментами отбора проб по мере изменения состава реагентов. Следовательно, скорость отвода тепла превышает указанный максимум 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 и Time-variable 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] Шринивасан, Б., С. Паланки и Д. Бонвин, «Динамическая оптимизация периодических процессов I. Характеристика номинального решения», Computers and Chemical Engineering, vol. 27 (2003), pp. 1-26.