Дистилляционная башня

Дистилляционная колонна использует 29 идеальных стадий для разделения смеси бензола, толуола и ксилолов (представленных п-ксилолом). Процесс дистилляции непрерывен. Оборудование включает ребойлер и полный конденсатор, как показано ниже:

Дистилляционная башня работает в номинальном установившемся условии:

Цели управления перечислены ниже, отсортированы по их важности:

Создайте высокоточную модель объекта управления в динамике Aspen Plus

Использование приложения Aspen Plus RADFRAC блок для определения статических характеристик башни. В дополнение к обычной информации, необходимой для статической симуляции, необходимо задать гидравлику поддона, геометрию отстойника башни и размер емкости орошения. Лотки представляют собой проект, разнесенную на 18 дюймов друг от друга. Все лотки имеют диаметр 1,95 м с высотой вейра 5 см. Номинальные глубины жидкости составляют 0,67 м и 1,4875 м в горизонтальной емкости орошения и отстойнике соответственно.

Установившаяся модель портируется на Aspen Plus Dynamics (APD) для управляемой потоком симуляции. Это пренебрегает динамикой привода и принимает точное регулирование управляемых скоростей потока жидкости. По умолчанию APD добавляет ПИ-контроллеров, чтобы регулировать давление в башне и два уровня жидкости. В этом примере ПИ-контроллеров по умолчанию намеренно удаляются.

Модель APD высокоточной дистилляционной башни показана ниже:

Линеаризация объекта с помощью интерфейса Aspen Plus Система Управления

Model Predictive Controller требует модель LTI объекта. В этом примере входы объекта управления:

Выходные выходы объекта:

Aspen Plus Dynamics предоставляет инструмент Система Управления Interface (CDI), который линеаризирует динамическую модель в заданном условии.

Следующие шаги сделаны, чтобы получить линейную модель объекта управления в Aspen Plus Dynamics.

Step 1: Добавьте скрипт к модели APD под Flowsheet папка. В этом примере имя скрипта CDI_Calcs (как показано выше) и содержит следующие команды APD:

Set Doc = ActiveDocument
set CDI = Doc.CDI
CDI.Reset
CDI.AddInputVariable "blocks(""B1"").condenser(1).QR"
CDI.AddInputVariable "blocks(""B1"").QrebR"
CDI.AddInputVariable "blocks(""B1"").Reflux.FmR"
CDI.AddInputVariable "streams(""2"").FmR"
CDI.AddInputVariable "streams(""3"").FmR"
CDI.AddInputVariable "streams(""1"").FR"
CDI.AddOutputVariable "blocks(""B1"").Stage(1).P"
CDI.AddOutputVariable "blocks(""B1"").Stage(1).Level"
CDI.AddOutputVariable "blocks(""B1"").SumpLevel"
CDI.AddOutputVariable "streams(""2"").Zmn(""TOLUENE"")"
CDI.AddOutputVariable "streams(""3"").Zmn(""BENZENE"")"
CDI.Calculate

Step 2: Инициализируйте модель APD до номинального установившегося условия.

Step 3: Активируйте скрипт, который генерирует следующие текстовые файлы:

В этом примере cdi_list.lis включает следующую информацию:

A matrix computed, number of non-zero elements = 1408
B matrix computed, number of non-zero elements = 26
C matrix computed, number of non-zero elements = 20
G matrix computed, number of non-zero elements = 30
Number of state variables:    120
Number of input variables:      6
Number of output variables:     5
Input variables:
   1        -3690034.247458334  BLOCKS("B1").Condenser(1).QR
   2            3819023.193875  BLOCKS("B1").QRebR
   3            22135.96620144  BLOCKS("B1").Reflux.FmR
   4            11717.39655353  STREAMS("2").FmR
   5            34352.86345834  STREAMS("3").FmR
   6                       500  STREAMS("1").FR
Output variables:
   1         1.100022977953499  BLOCKS("B1").Stage(1).P
   2        0.6700005140605662  BLOCKS("B1").Stage(1).Level
   3                    1.4875  BLOCKS("B1").SumpLevel
   4       0.05002582161855798  STREAMS("2").Zmn("TOLUENE")
   5       0.01705308738356429  STREAMS("3").Zmn("BENZENE")

Номинальные значения переменных состояния, перечисленных в файле, игнорируются, поскольку они не требуются в проекте MPC.

Создайте масштабированную и уменьшенную модель пространства состояний LTI

Step 1: Преобразуйте CDI-сгенерированные разреженные матрицы в модель пространства состояний.

Загрузите матрицы пространства состояний из файлов данных CDI в рабочее пространство MATLAB и преобразуйте разреженные матрицы в полные матрицы.

load mpcdistillation_cdi_A.dat
load mpcdistillation_cdi_B.dat
load mpcdistillation_cdi_C.dat
A = full(spconvert(mpcdistillation_cdi_A));
B = full(spconvert(mpcdistillation_cdi_B));
C = full(spconvert(mpcdistillation_cdi_C));
D = zeros(5,6);

Возможно, что вся разреженная матричная строка или столбец равны нулю, и в этом случае вышеуказанных команд недостаточно. Используйте следующие дополнительные проверки, чтобы убедиться, что A, B и C имеют правильные размерности:

[nxAr,nxAc] = size(A);
[nxB,nu] = size(B);
[ny,nxC] = size(C);
nx = max([nxAr, nxAc, nxB, nxC]);
if nx > nxC
    C = [C, zeros(ny,nx-nxC)];
end
if nx > nxAc
    A = [A zeros(nxAr,nx-nxAc)];
end
if nx > nxAr
    nxAc = size(A,2);
    A = [A; zeros(nx-nxAr, nxAc)];
end
if nxB < nx
    B = [B; zeros(nx-nxB,nu)];
end

Step 2: Масштабирование сигналов объекта.

Хорошая практика, если не существенная, преобразовать сигналы объекта от технических модулей в однородную безразмерную шкалу (например, 0-1 или 0-100%). Одной из альтернативных вариантов является определение масштабных коэффициентов как части проекта Model Predictive Controller. Это может значительно упростить настройку контроллера. Для получения примера смотрите, Использование Шкалы факторов для облегчения настройки веса.

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

U_span = [2*(-3690034), 2*3819023, 2*22136, 2*11717, 2*34353, 2*500];
U_nom = 0.5*U_span;
U_zero = zeros(1,6);
Y_nom = [1.1, 0.67, 1.4875, 0.050026, 0.017053];
Y_span = [0.4, 2*Y_nom(2:5)];
Y_zero = [0.9, 0, 0, 0, 0];

Масштабируйте матрицы B и C так, чтобы все входные/выходные переменные были выражены в процентах.

B = B.*(ones(nx,1)*U_span);
C = C./(ones(nx,1)*Y_span)';

Step 3: Определите модель объекта управления в пространстве состояний.

G = ss(A,B,C,D);
G.TimeUnit = 'hours';
G.u = {'Qc','Qr','R','D','B','F'};
G.y = {'P','RLev','Slev','xD','xB'};

Step 4: Уменьшите порядок модели.

Снижение сложности модели ускоряет вычисления с незначительным эффектом на точность предсказания. Используйте hsvd команда для определения состояний, которые могут быть безопасно сброшены. Используйте balred функция, чтобы удалить эти состояния и уменьшить порядок модели.

[hsv, baldata] = hsvd(G);
order = find(hsv>0.01,1,'last');
Options = balredOptions('StateElimMethod','Truncate');
G = balred(G,order,baldata,Options);

Исходная модель имеет 120 состояний, а уменьшенная модель имеет только 16 состояний. Обратите внимание, что Truncate опция используется в balred функция для сохранения нуля D- матрицы. Модель имеет два полюса в нуле, которые соответствуют двум уровням жидкости.

Точность тестирования линейной модели объекта управления

Перед продолжением проекта MPC рекомендуется проверить, что масштабированная модель LTI точна для небольших изменений в входах объекта. Для этого необходимо сравнить ответ нелинейного объекта в APD и ответ линейной модели G.

Step 1: Чтобы получить ответ нелинейного объекта, создайте модель Simulink и добавьте к ней блок Aspen Modeler.

Блок предоставляется Aspen Plus Dynamics в своих AMSimulink библиотека.

Step 2: Дважды кликните блок и укажите местоположение модели APD.

Информация о модели APD затем импортируется в Simulink. Для больших моделей APD процесс импорта может занять некоторое время.

Step 3: Задайте входной и выходной сигналы в блоке AMSimulation.

Используйте те же имена переменных и последовательность, что и в скрипте CDI.

Теперь блок показывает входные и выходные порты для каждого заданного сигнала.

Step 4: Разверните модель Simulink с входным сигналом, поступающим от переменной Umat и выходной сигнал, сохраненный в переменной Ypct_NL. Обе переменные создаются в шаге 5.

Начиная с Umat находится в процентных модулях, Pct2Engr блок реализован для преобразования из процентных модулей в инженерные модули.

Начиная с Ypct_NL в процентных модулях реализован блок «Engr2Pct» для преобразования из технических модулей в процентные модули.

Со всем подключенным и настроенным модель выглядит следующим образом:

Step 5: Проверьте линейную модель с косимуляцией.

В этом примере в качестве сигнала возбуждения для объекта используется увеличение масштабируемой скорости флегмы на 1 процент (вход # 3 ).

U_nom_pct = (U_nom - U_zero)*100./U_span;   % Convert nominal condition from engineering units to percentages
Y_nom_pct = (Y_nom - Y_zero)*100./Y_span;
Tend = 1;                                   % Simulation duration (1 hour)
t = (0:1/60:Tend)';                         % Sample period is 1 minute
nT = length(t);
Upct = ones(nT,1)*U_nom_pct;
DUpct = zeros(nT,6);
DUpct(:,3) = ones(nT,1);                    % Input signal where step occurs in channel #3

Реакция линейной модели объекта управления вычисляется с помощью lsim команда и сохранена в переменных Ypct_L.

Ypct_L = lsim(G,DUpct,t);
Ypct_L = Ypct_L + ones(nT,1)*Y_nom_pct;

Реакция нелинейного объекта получена путем косимуляции между Simulink и Aspen Plus Dynamics. Сигнал возбуждения Umat построено следующим образом. Результат хранится в переменной Ypct_NL.

Umat = [t, Upct+DUpct];

Сравните линейные и нелинейные отклики модели.

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

Можно повторить вышеописанные шаги, чтобы проверить аналогичное согласие для других пяти входов.

Модель проекта

Учитывая модель предсказания LTI, вы готовы к разработке прогнозирующего контроллера модели. В этом примере манипулируемые переменные являются первыми пятью входами объекта. Шестой вход объекта (feed скорости потока жидкости) является измеренным нарушением порядка для feedforward компенсации. Все выходы объекта измеряются.

Step 1: Увеличить объект, чтобы смоделировать неизмеренные нарушения порядка нагрузки.

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

В этом примере 5 неизмеренных нарушений порядка нагрузки добавляются к модели объекта управления G. Всего в настоящее время существует 11 входов в модель предсказания Gmpc: 5 манипулированных переменных, 1 измеренное нарушение порядка и 5 неизмеренных нарушения порядка.

Gmpc = ss(G.A,G.B(:,[1:6,1:5]),G.C,zeros(5,11),'TimeUnit','hours');
InputName = cell(1,11);
for i = 1:5
    InputName{i} = G.InputName{i};
    InputName{i+6} = [G.InputName{i}, '-UD'];
end
InputName{6} = G.InputName{6};
Gmpc.InputName = InputName;
Gmpc.InputGroup = struct('MV',1:5,'MD',6,'UD',7:11);
Gmpc.OutputName = G.OutputName;

Step 2: Создайте начальный прогнозирующий контроллер модели и задайте шаг расчета и горизонты.

В этом примере период дискретизации контроллера составляет 30 секунд. Горизонт предсказания является 60 интервалами (30 минуты), что достаточно велико, чтобы сделать эффективность контроллера нечувствительным к дальнейшим увеличениям горизонта предсказания. Горизонт управления является 4 интервалами (2 минуты), что относительно мало для уменьшения вычислительных усилий.

Ts = 30/3600;       % sample time
PH = 60;            % prediction horizon
CH = 4;             % control horizon 
MPCobj = mpc(Gmpc,Ts,PH,CH);  % MPC object

-->The "Weights.ManipulatedVariables" property of "mpc" object is empty. Assuming default 0.00000.
-->The "Weights.ManipulatedVariablesRate" property of "mpc" object is empty. Assuming default 0.10000.
-->The "Weights.OutputVariables" property of "mpc" object is empty. Assuming default 1.00000.

Step 3: Задайте веса для управляемых переменных и управляемых выходов.

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

Нет причин держать конкретное СН в заданном значении, поэтому установите Weights.ManipulatedVariables свойство в нуль:

MPCobj.Weights.ManipulatedVariables = [0, 0, 0, 0, 0];

Продукт дистиллята (MV # 4) поступает на хранение. Единственным СН, влияющим на работу установки ниже по потоку, является скорость забоя (СН # 5). Чтобы предотвратить быстрые изменения скорости дна, сохраните вес по умолчанию 0,1 для его скорости изменения. Уменьшите другую скорость изменения веса в 10 раз:

MPCobj.Weights.ManipulatedVariablesRate = [0.01, 0.01, 0.01, 0.01, 0.1];

Цели управления обеспечивают рекомендации по выбору весов на контролируемых выходах:

Учитывая эти приоритеты, веса на контролируемых выходах выбираются следующим образом:

MPCobj.Weights.OutputVariables = [10, 0.1, 0.1, 1, 0.5];

Масштабирование модели упрощает выбор оптимизационных весов. В противном случае, в дополнение к относительному приоритету каждой переменной, вы также должны будете рассмотреть относительные величины переменных и выбрать веса соответственно.

Step 4: Задайте номинальные входные/выходные значения объекта управления.

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

MPCobj.Model.Nominal.U = [U_nom_pct'; zeros(5,1)];
MPCobj.Model.Nominal.Y = Y_nom_pct';

Step 5: Настройте коэффициент усиления оценщика состояния.

Регулировка усиления оценщика состояния влияет на эффективность подавления помех. Увеличение усиления оценщика состояния (например, путем увеличения усиления модели входного/выходного возмущения) заставляет контроллер более агрессивно реагировать на выходные изменения (потому что контроллер принимает, что основным источником выходных изменений является нарушение порядка, а не шум измерения). С другой стороны, уменьшение усиления устройства оценки состояния делает систему с обратной связью более устойчивой.

Во-первых, проверьте, обеспечивает ли использование оценки состояния по умолчанию достойное подавление помех эффективности.

Симулируйте реакцию с обратной связью на шаг 1% модулей в орошении (MV # 3) в MATLAB. В симуляции используются G как объект, что не подразумевает несоответствия модели.

T = 30;                                 % Simulation time
r = Y_nom_pct;                          % Nominal setpoints
v = U_nom_pct(6);                       % No measured disturbance
SimOptions = mpcsimopt(MPCobj);
SimOptions.InputNoise = [0 0 1 0 0];    % 1% unit step in reflux
[y_L,t_L,u_L] = sim(MPCobj, T, r, v, SimOptions); % Closed-loop simulation

-->Converting model to discrete time.
-->The "Model.Disturbance" property of "mpc" object is empty:
   Assuming unmeasured input disturbance #7 is integrated white noise.
   Assuming unmeasured input disturbance #8 is integrated white noise.
   Assuming unmeasured input disturbance #9 is integrated white noise.
   Assuming unmeasured input disturbance #10 is integrated white noise.
   Assuming unmeasured input disturbance #11 is integrated white noise.
   Assuming no disturbance added to measured output channel #1.
   Assuming no disturbance added to measured output channel #4.
   Assuming no disturbance added to measured output channel #5.
   Assuming no disturbance added to measured output channel #2.
   Assuming no disturbance added to measured output channel #3.
-->The "Model.Noise" property of the "mpc" object is empty. Assuming white noise on each measured output channel.

% plot responses
f1 = figure();
subplot(2,1,1);
plot(t_L,y_L,[0 t_L(end)],[50 50],'k--')
title('Controlled Outputs, %')
legend(Gmpc.OutputName,'Location','NorthEastOutside')
subplot(2,1,2);
plot(t_L,u_L(:,1:5),[0 t_L(end)],[50 50],'k--')
title('Manipulated Variables, %')
legend(Gmpc.InputName(1:5),'Location','NorthEastOutside')
xlabel('Time, h')

Оценка по умолчанию обеспечивает медленное отклонение нагрузки. В частности, критические xD вывод падает до 49% и только начал возвращаться к уставке через 0,25 часов.

Во-вторых, увеличьте коэффициент усиления оценщика путем умножения коэффициента усиления модели входного возмущения по умолчанию в 25 раз.

EstGain = 25;                       % factor of 25
Gd = getindist(MPCobj);             % get default input disturbance model
Gd_new = EstGain*Gd;                % create new input disturbance model
setindist(MPCobj,'Model',Gd_new); % set input disturbance model
[y_L,t_L,u_L] = sim(MPCobj,T,r,v,SimOptions); % Closed-loop simulation

-->Converting model to discrete time.
   Assuming no disturbance added to measured output channel #1.
   Assuming no disturbance added to measured output channel #4.
   Assuming no disturbance added to measured output channel #5.
   Assuming no disturbance added to measured output channel #2.
   Assuming no disturbance added to measured output channel #3.
-->The "Model.Noise" property of the "mpc" object is empty. Assuming white noise on each measured output channel.

% plot responses
f2 = figure();
subplot(2,1,1);
plot(t_L,y_L,[0 t_L(end)],[50 50],'k--')
title('Controlled Outputs, %')
legend(Gmpc.OutputName,'Location','NorthEastOutside')
subplot(2,1,2)
plot(t_L,u_L(:,1:5),[0 t_L(end)],[50 50],'k--')
title('Manipulated Variables, %')
legend(Gmpc.InputName(1:5),'Location','NorthEastOutside')
xlabel('Time, h')

Теперь пиковое отклонение в xD на 50% меньше, чем в случае по умолчанию и xD возвращается к своей уставке намного быстрее. Другие переменные также реагируют быстрее.

В-третьих, посмотрите на рефлюксную реакцию (# 3 на графике «Манипулированные переменные»). Поскольку нарушение порядка является шагом в 1% модулей, реакция начинается с 51%, и его окончательное значение составляет 50% в установившемся состоянии. Рефлюксная реакция перерегулировала на 20% (достигая 49,8%) перед осаждением. Такое количество перерегулирования приемлемо.

Если бы коэффициент усиления оценщика был увеличен дополнительно (например, в 50 раз), коэффициент перерегулирования контроллера тоже увеличился бы. Однако такое агрессивное поведение вряд ли будет устойчивым при применении к нелинейной модели объекта управления.

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

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

Как правило, далее проверяется реакция на изменения уставки. Если реакция слишком агрессивна, можно использовать фильтр уставки, чтобы сглаживать его. Фильтр уставки не влияет на подавление помех нагрузки и, таким образом, может быть настроен независимо.

Косимуляция контроллера MPC и нелинейного объекта

Используйте косимуляцию, чтобы определить, является ли проект MPC достаточно устойчивым, чтобы управлять нелинейной моделью объекта управления.

Step 1: Добавьте ограничения к контроллеру MPC

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

MV = MPCobj.MV;
OV = MPCobj.OV;
% Physical bounds on MVs at 0 and 100 
for i = 1:5
    MV(i).Min = 0;
    MV(i).Max = 100;
end
MPCobj.MV = MV;
% Keep liquid levels greater than 25% and less than 75% of capacity.
for i = 2:3
    OV(i).Min = 25;
    OV(i).Max = 75;
end
MPCobj.OV = OV;

Step 2: Создайте модель Simulink для косимуляции.

Модель может моделировать шаг 1% модулей в рефлюксе (MV # 3). Это также может имитировать изменение состава сырья, которое является общим нарушением порядка и отличается от нарушений порядка нагрузки, явно рассматриваемых в проекте.

Step 3: Моделируйте шаг 1% модулей в орошении (MV # 3). Сравните отклики замкнутой системы между использованием линейной модели объекта управления и нелинейной модели объекта управления.

График состава продукта дистиллята (xD) и скорость рефлюкса (R):

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

Step 4: Моделируйте большое уменьшение фракции бензола (от 0,3 до 0,22) в потоке сырья. Сравните отклики с обратной связью между использованием линейной и нелинейной моделей объекта управления.

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