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

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

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

Цели контроля перечислены ниже с разбивкой по их важности:

Построение модели завода высокой четкости в Aspen Plus Dynamics

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

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

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

Линеаризация установки с использованием интерфейса проектирования Aspen Plus Control

Модельный предиктивный контроллер требует модели завода LTI. В этом примере входные данные установки:

Выходы установки:

Aspen Plus Dynamics предоставляет инструмент интерфейса проектирования управления (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%). Одним из альтернативных вариантов является определение коэффициентов масштаба как части конструкции контроллера прогнозирования модели. Это значительно упрощает настройку контроллера. Например, см. раздел Использование масштабных коэффициентов для облегчения настройки веса.

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

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% масштабированной скорости орошения (вход No 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, вы готовы разработать контроллер прогнозирования модели. В этом примере управляемые переменные являются первыми пятью входными данными установки. Шестой вход установки (расход сырья) представляет собой измеренное возмущение для компенсации потока сырья. Измеряются все выходы установки.

Step 1Увеличение мощности установки для моделирования неизмеренных нарушений нагрузки.

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

В этом примере к модели завода добавляется 5 неизмеренных нарушений нагрузки. G. В общей сложности сейчас существует 11 входных данных для модели прогнозирования Gmpc5 манипулируемых переменных, 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) поступает на хранение. Единственным MV, влияющим на работу установки в нисходящем направлении, является кубовая скорость (MV # 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Укажите номинальные значения ввода/вывода для завода.

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

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 раз), перерасход контроллера также увеличился бы. Однако такое агрессивное поведение вряд ли будет устойчивым при применении к нелинейной модели установки.

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

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

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

Контроллер Cosimulate MPC и нелинейная установка

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

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

Поскольку нелинейная модель установки имеет ограничения на вход и выход во время работы, ограничения на СН и 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 часов.