Exponenta Banner
exponenta event banner

Проектное исследование с использованием выборки параметров (код)

Открыть сценарий

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

Анализ конструкции CSTR осуществляется путем определения параметров конструкции с использованием распределения вероятностей. Распределения используются для генерации случайных выборок в пространстве проектирования и выполнения оценки конструкции по методу Монте-Карло в этих точках выборки. Затем создаются графики для визуализации пространства проектирования и выбирается наилучшая конструкция. Затем можно использовать наилучшую конструкцию в качестве первоначального предположения для оптимизации конструкции.

Можно также использовать выборку расчетного пространства и результаты оценки Монте-Карло для анализа влияния проектных параметров на проект, см. раздел Определение ключевых параметров для оценки (код).

Модель реактора емкости с непрерывным перемешиванием (CSTR)

Реакторы резервуаров с непрерывным перемешиванием (CSTR) являются обычными в технологической промышленности. модель Simulink, sdoCSTR, моделирует реактор с диабетической оболочкой (т.е. неадиабатический), описанный в [1]. Предполагается, что CSTR идеально смешивается с одной экзотермической и необратимой реакцией первого порядка., $A\rightarrow B$$A$реагент, превращается $B$в продукт.

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

$$\frac{d C_A}{dt} = \frac{F}{A * h}(C_{feed} - C_A)-r*C_A$$

$$\frac{d T}{dt} = \frac{F}{A * h}(T_{feed} - T) - \frac{H}{c_p \rho} r
-\frac{U}{c_p*\rho*h}(T - T_{cool})$$

$$r = k_0*e^{\frac{-E}{R*T}}$$

  • $C_A$и$C_{feed}$ - Концентрации А в CSTR и в сырье [кгмоль/м ^ 3]

  • $T$, $T_{feed}$и -$T_{cool}$ CSTR, температура сырья и хладагента [K]

  • $F$ и$\rho$ - Объемный расход [м ^ 3/ч] и плотность материала в CSTR [1/м ^ 3]

  • $h$ и$A$ - Высота [м] и площадь нагретого поперечного сечения [м ^ 2] CSTR.

  • $k_0$ - Предэкспоненциальный нетермический коэффициент для реакции$A\rightarrow B$ [1/ч]

  • $E$ и$H$ - Энергия активации и теплота реакции для [$A\rightarrow B$ккал/кгмоль]

  • $R$ - газовая постоянная Больцмана [ккал/( кгмоль * К)]

  • $c_p$ и$U$ - теплоемкость [ккал/К] и коэффициенты теплопередачи [ккал/( м ^ 2 * К * ч)]

Откройте модель Simulink.

open_system('sdoCSTR');

Проблема проектирования CSTR

Предположим, что CSTR имеет цилиндрическую форму, а хладагент приложен к основанию цилиндра. Настройте площадь поперечного сечения CSTR $A$и высоту CSTR $h$для достижения следующих целей проектирования:

  • Минимизируйте изменение остаточной концентрации,. $C_A$Колебания остаточной концентрации негативно влияют на качество продукта CSTR. Минимизация изменений также повышает прибыль CSTR.

  • Минимизируйте среднюю температуру хладагента. $T_{cool}$Нагрев или охлаждение охлаждающей жидкости рубашки обходится дорого. Минимизация средней температуры хладагента повышает прибыль от CSTR.

Конструкция должна обеспечивать возможность изменения качества концентрации подаваемого сырья и $C_{feed}$температуры сырья. $T_{feed}$В CSTR подается сырье от разных поставщиков. Качество сырья отличается от поставщика к поставщику, а также варьируется в пределах каждой партии поставки.

Задание конструктивных переменных

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

  • Площадь поперечного сечения цилиндра $A$

  • Высота цилиндра $h$

p = sdo.getParameterFromModel('sdoCSTR',{'A','h'});

Ограничьте площадь поперечного сечения диапазоном [0,2 2] м ^ 2.

p(1).Minimum = 0.2;
p(1).Maximum = 2;

Ограничьте высоту диапазоном [0,5 3] м.

p(2).Minimum = 0.5;
p(2).Maximum = 3;

Образец пространства проектирования

Создайте пространство параметров для конструктивных переменных. Пространство параметров характеризует допустимые значения параметров и комбинации значений параметров.

pSpace = sdo.ParameterSpace(p);

В пространстве параметров для конструктивных переменных используются равномерные распределения по умолчанию. Нижние и верхние границы распределения задаются минимальными и максимальными значениями конструктивных переменных соответственно.

Используйте sdo.sample для создания выборок из пространства параметров. Образцы используются для оценки модели и исследования пространства проектирования.

rng default; % For reproducibility
pSmpl = sdo.sample(pSpace,30);

Используйте sdo.scatterPlot для визуализации пространства параметров выборки. График рассеяния показывает распределения параметров на диагональных субплотах и попарные комбинации параметров на внедиагональных субплотах.

figure, sdo.scatterPlot(pSmpl)

Указать неопределенные переменные

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

pUnc = sdo.getParameterFromModel('sdoCSTR',{'FeedCon0','FeedTemp0'});

Создайте пространство параметров для неопределенных переменных. Используйте нормальные распределения для обеих переменных. Укажите среднее значение в качестве текущего значения параметра. Укажите дисперсию 5% от среднего значения для концентрации сырья и 1% от среднего значения для температуры.

uSpace = sdo.ParameterSpace(pUnc);
uSpace = setDistribution(uSpace,'FeedCon0',makedist('normal',pUnc(1).Value,0.05*pUnc(1).Value));
uSpace = setDistribution(uSpace,'FeedTemp0',makedist('normal',pUnc(2).Value,0.01*pUnc(2).Value));

Концентрация сырья обратно коррелирует с температурой сырья. Добавьте эту информацию в пространство параметров.

uSpace.RankCorrelation = [1 -0.6; -0.6 1];

Матрица ранговой корреляции имеет строку и столбец для каждого параметра с элементом (i, j), определяющим корреляцию между параметрами i и j.

Образец пространства параметров. График рассеяния показывает корреляцию между концентрацией и температурой.

uSmpl = sdo.sample(uSpace,60);
sdo.scatterPlot(uSmpl)

В идеале нужно оценить конструкцию для каждой комбинации точек в конструкции и неопределенных пространств, что подразумевает 30 * 60 = 1800 моделирования. Каждое моделирование занимает около 0,5 сек. Вы можете использовать параллельные вычисления для ускорения оценки. В этом примере вместо этого используются только образцы с максимальными и минимальными значениями концентрации и температуры, что сокращает время оценки примерно до 1 мин.

[~,iminC] = min(uSmpl.FeedCon0);
[~,imaxC] = max(uSmpl.FeedCon0);
[~,iminT] = min(uSmpl.FeedTemp0);
[~,imaxT] = max(uSmpl.FeedTemp0);
uSmpl = uSmpl(unique([iminC,imaxC,iminT,imaxT]) ,:)

uSmpl =

  4x2 table

    FeedCon0    FeedTemp0
    ________    _________

     9.4555      303.58  
     11.175      288.13  
     11.293      290.73  
     8.9308      294.16

Создание функции оценки

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

Определение требований к конструкции

Оценка точки в пространстве проектирования требует регистрации сигналов модели. Регистрируемые сигналы используются для оценки проектных требований.

Запишите следующие сигналы:

  • Концентрация CSTR, доступная на втором выходном порту sdoCSTR/CSTR блок

Conc = Simulink.SimulationData.SignalLoggingInfo;
Conc.BlockPath               = 'sdoCSTR/CSTR';
Conc.OutputPortIndex         = 2;
Conc.LoggingInfo.NameMode    = 1;
Conc.LoggingInfo.LoggingName = 'Concentration';

  • Температура теплоносителя, доступная на первом выходе sdoCSTR/Controller блок

Coolant = Simulink.SimulationData.SignalLoggingInfo;
Coolant.BlockPath               = 'sdoCSTR/Controller';
Coolant.OutputPortIndex         = 1;
Coolant.LoggingInfo.NameMode    = 1;
Coolant.LoggingInfo.LoggingName = 'Coolant';

Создайте и настройте имитационный тестовый объект для регистрации требуемых сигналов.

simulator = sdo.SimulationTest('sdoCSTR');
simulator.LoggingInfo.Signals = [Conc,Coolant];

Функция оценки

Используйте анонимную функцию с одним аргументом, вызывающим sdoCSTR_design функция.

evalDesign = @(p) sdoCSTR_design(p,simulator,pUnc,uSmpl);

evalDesign функция:

  • Имеет один входной аргумент, указывающий размеры CSTR

  • Возвращает целевое значение оптимизации

sdoCSTR_design функция использует for цикл, который итерирует значения пробы, указанные для концентрации и температуры сырья. В цикле функция:

  • Моделирование модели с использованием текущей расчетной точки, концентрации сырья и значений температуры сырья

  • Расчет изменения остаточной концентрации и затрат на температуру хладагента

Для просмотра целевой функции введите edit sdoCSTR_design.

type sdoCSTR_design

function design = sdoCSTR_design(p,simulator,pUnc,smplUnc)
%SDOCSTR_DESIGN
%
% The sdoCSTR_design function is used to evaluate a CSTR  design.
%
% The |p| input argument is the vector of CSTR dimensions.
%
% The |simulator| input argument is a sdo.SimulinkTest object used to
% simulate the |sdoCSTR| model and log simulation signals.
%
% The |pUnc| input argument is a vector of parameters to specify the CSTR
% input feed concentration and feed temperature. The |smplUnc| argument is
% a table of different feed concentration and temperature values.
%
% The |design| return argument contains information about the design
% evaluation that can be used by the |sdo.optimize| function to optimize
% the design.
%
% see also sdo.optimize, sdoExampleCostFunction
%

% Copyright 2012-2013 The MathWorks, Inc.

%% Model Simulations and Evaluations
%
% For each value in |smplUnc|, configure and simulate the model. Use
% the logged concentration and coolant signals to compute the design cost.
%
costConc    = 0;
costCoolant = 0;
for ct=1:size(smplUnc,1)
    %Set the feed concentration and temperature values
    pUnc(1).Value = smplUnc{ct,1};
    pUnc(2).Value = smplUnc{ct,2};
    
    %Simulate model
    simulator.Parameters = [p; pUnc];
    simulator = sim(simulator);
    logName   = get_param('sdoCSTR','SignalLoggingName');
    simLog    = get(simulator.LoggedData,logName);
    
    %Compute Concentration cost based on the standard deviation of the
    %concentration from a nominal value.
    Sig = find(simLog,'Concentration');
    costConc = costConc+10*std(Sig.Values-2);
    
    %Compute coolant cost based on the mean deviation from room
    %temperature.
    Sig = find(simLog,'Coolant');
    costCoolant = costCoolant+abs(mean(Sig.Values - 294))/30;
end

%% Return Total Cost 
%
% Compute the total cost as a sum of the concentration and coolant costs.
%
design.F = costConc + costCoolant;

%%
% Add the individual cost terms to the return argument. These are not used
% by the optimizer, but included for convenience.
design.costConc    = costConc;               
design.costCoolant = costCoolant;
end

Оценить

Используйте sdo.evaluate для оценки модели в точках проектирования образца.

y = sdo.evaluate(evalDesign,p,pSmpl);

Model evaluated at 30 samples.

Просмотрите результаты оценки с помощью графика разброса. График рассеяния показывает попарные графики для каждой конструктивной переменной (A, h) и проектной стоимости. Участок включает в себя общую стоимость ,F, а также затраты на теплоноситель и концентрацию, costCoolant и costConc соответственно.

sdo.scatterPlot(pSmpl,y);

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

Создайте сетчатый график, показывающий общую стоимость как функцию A и h.

Ftotal = scatteredInterpolant(pSmpl.A,pSmpl.h,y.F);
xR = linspace(min(pSmpl.A),max(pSmpl.A),60);
yR = linspace(min(pSmpl.h),max(pSmpl.h),60);
[xx,yy] = meshgrid(xR,yR);
zz = Ftotal(xx,yy);
mesh(xx,yy,zz)
view(56,30)
title('Total cost as function of A and h')
zlabel('Ftotal')
xlabel(p(1).Name), ylabel(p(2).Name);

График показывает, что высокие значения A и h приводят к снижению затрат. Наилучшая конструкция в отобранном пространстве соответствует конструкции с наименьшей стоимостью.

[~,idx] = min(y.F);
pBest = [y(idx,:), pSmpl(idx,:)]

pBest =

  1x5 table

      F      costConc    costCoolant      A         h   
    _____    ________    ___________    ______    ______

    2.106     1.5505       0.55552      1.9271    2.3867

Уточнение пространства проектирования

График общей стоимости поверхности показывает, что проекты с низкой стоимостью представляют собой проекты с A в диапазоне [1.5 2] и h в диапазоне [2 3]. Измените распределения пространства параметров для A и h и выполните повторную выборку пространства проектирования, чтобы сосредоточиться на этой области.

pSpace = setDistribution(pSpace,'A',makedist('uniform',1.5,2));
pSpace = setDistribution(pSpace,'h',makedist('uniform',2,3));
pSmpl = sdo.sample(pSpace,30);

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

pSmpl = [pSmpl;pBest(:,4:5)];
sdo.scatterPlot(pSmpl)

Оценка с использованием уточненного пространства проектирования

y = sdo.evaluate(evalDesign,p,pSmpl);

Model evaluated at 31 samples.

Создайте сеточный график для этого раздела пространства проектирования. Поверхность показывает, что лучшие конструкции находятся вблизи точки A = 1,9, h = 2,1.

Ftotal = scatteredInterpolant(pSmpl.A,pSmpl.h,y.F);
xR = linspace(min(pSmpl.A),max(pSmpl.A),60);
yR = linspace(min(pSmpl.h),max(pSmpl.h),60);
[xx,yy] = meshgrid(xR,yR);
zz = Ftotal(xx,yy);
mesh(xx,yy,zz)
view(56,30)
title('Total cost as function of A and h')
zlabel('Ftotal')
xlabel(p(1).Name), ylabel(p(2).Name);

Найдите лучший проект из нового пространства проектирования и сравните его с лучшей точкой проектирования, найденной ранее.

[~,idx] = min(y.F);
pBest = [pBest; [y(idx,:), pSmpl(idx,:)]]

pBest =

  2x5 table

      F       costConc    costCoolant      A         h   
    ______    ________    ___________    ______    ______

     2.106     1.5505       0.55552      1.9271    2.3867
    1.9754     1.4824       0.49295      1.9695    2.1174

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

Связанные примеры

Сведения о том, как исследовать пространство проектирования CSTR с помощью анализатора чувствительности, см. в разделе Исследование конструкции с использованием выборки параметров (GUI).

Как оптимизировать конструкцию CSTR с помощью sdo.optimize см. раздел Оптимизация конструкции с неопределенными переменными (код).

Изучение способов анализа влияния параметров конструкции на конструкцию с помощью sdo.analyze , см. раздел Определение ключевых параметров для оценки (код)

Ссылки

[1] Becette, B.W. Process Dynamics: моделирование, анализ и моделирование. 1-я ред. река Верхнее Седло, Нью-Джерси: Прентис Холл, 1998.

Закрыть модель

bdclose('sdoCSTR')

Связанные темы

Документация по оптимизации проектирования Simulink

Поддержка