Документация

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

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

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

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

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

Переменные, которыми управляют,

Измеренное воздействие

Состояния

Выходные параметры

Цель состоит в том, чтобы максимизировать производство 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.

Анализ результатов Optmization

Найдите оптимальные траектории для переменных, которыми управляют, таким образом, что производство 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 = 6.56014 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.55228 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.31788 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.