Пакетный ФРС химический реактор

Следующие необратимые и экзотермические реакции происходят в пакетном реакторе [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)')

Тщательное изучение уровня удаления тепла показывает, что может показать 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 = 3.20304 sec. Production of C = 2.01754 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.99957 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.54985 sec. Production of C = 2.02349 mol

В большинстве приложений MPC линейные решения MPC, такие как Адаптивный 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] Сринивасан, B., С. Паланки и Д. Бонвин, "Динамическая оптимизация процессов пакетной обработки I. Характеристика номинального решения", Компьютеры и Химическое машиностроение, издание 27 (2003), стр 1-26.