Модель Simulink

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

mdl = 'scdTanks';
open_system(mdl)

Желателен поток между Tank1 и Tank2. Поток между Tank2 и Tank3 является нежелательной неизбежной утечкой.

При ожидаемом устойчивом состоянии этой системы:

Из-за слабой связности между Tank1 и Tank2 трудно обрезать модель так, чтобы давления в Tank1 и Tank2 были равными.

Обрезка модели без настроек

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

opspec = operspec(mdl);

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

opt = findopOptions('OptimizerType','lsqnonlin');
opt.DisplayReport = 'off';

Обрезать модель и просмотреть давления в усеченном резервуаре.

[op0,rpt0] = findop(mdl,opspec,opt);
op0.States

(1.) scdTanks/Inertia
      x: 0            
(2.) scdTanks/Tank1
      x: 9            
(3.) scdTanks/Tank2
      x: 9.5          
(4.) scdTanks/Tank3
      x: 10.5         

Урезанные давления в Tank1 и Tank2 не совпадают. Таким образом, спецификация рабочей точки по умолчанию не может найти рабочую точку, удовлетворяющую ожидаемым требованиям к установившемуся состоянию. При уменьшении допуска зависимости opt.OptimizationOptions.TolCon, невозможно достичь приемлемого стационарного решения из-за утечки между Tank2 и Tank3.

Добавить пользовательские ограничения

Чтобы задать пользовательские ограничения, определите функцию в текущей рабочей папке или на пути MATLAB с входными аргументами:

и выходные аргументы:

Каждый элемент c_ineq и c_eq задает одно ограничение. Определите определенные ограничения для приложения как алгебраические комбинации состояний, входов и выходов. Если нет пользовательских ограничений равенства или неравенства, верните соответствующий выходной аргумент как [].

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

function [c_ineq,c_eq] =  myConstraints(x,u,y)
    c_ineq = [];
    c_eq = [x(2)-x(3);     % Tank1 pressure - Tank2 pressure
            x(3)-x(4)+1];  % Tank2 pressure - Tank3 pressure + 1
end

Первая запись c_eq ограничивает давление Tank1 и Tank2 одним и тем же значением. Второе ограничение равенства определяет перепад давления между Tank2 и Tank3.

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

opspec.CustomConstrFcn = @myConstraints;

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

[op1,rpt1] = findop(mdl,opspec,opt);
op1.States

(1.) scdTanks/Inertia
      x: 0            
(2.) scdTanks/Tank1
      x: 9.33         
(3.) scdTanks/Tank2
      x: 9.33         
(4.) scdTanks/Tank3
      x: 10.3         

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

Чтобы проверить конечные значения указанных ограничений, можно проверить CustomEqualityConstr и CustomInequalityConstr свойства отчета о поиске в рабочей точке.

rpt1.CustomEqualityConstr

ans =

   1.0e-06 *

   -0.0001
   -0.1540

Значения, близкие к нулю, указывают, что ограничения равенства выполнены.

Добавить пользовательскую целевую функцию

Чтобы задать пользовательскую целевую функцию, определите функцию с теми же входными аргументами, что и пользовательская функция ограничения (x, u, и y) и выходной аргумент F. F - значение целевой функции, которое должно быть минимизировано во время обрезки, возвращаемое как скаляр.

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

В этом примере предположим, что давление в Tank3 должно находиться в диапазоне [16,20]. Однако это условие не всегда осуществимо. Поэтому, вместо того чтобы накладывать жесткие ограничения, добавьте объективную функцию, чтобы понести штраф, если давление не находится в [16,20] диапазоне. Для этого определите следующую пользовательскую целевую функцию.

function F =  myObjective(x,u,y)
    F = max(x(4)-20, 0) + max(16-x(4), 0);
end

Добавьте пользовательскую целевую функцию к объекту спецификации операционной точки.

opspec.CustomObjFcn = @myObjective;

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

[op2,rpt2] = findop(mdl,opspec,opt);
op2.States

(1.) scdTanks/Inertia
      x: 0            
(2.) scdTanks/Tank1
      x: 15           
(3.) scdTanks/Tank2
      x: 15           
(4.) scdTanks/Tank3
      x: 16           

В усеченной рабочей точке давление в Tank3 находится в диапазоне [16,20], указанном в пользовательской целевой функции.

Чтобы просмотреть конечное значение скалярной целевой функции, проверьте CustomObj свойства отчета о поиске в рабочей точке.

rpt2.CustomObj

ans =

     0

Добавить пользовательское сопоставление

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

Чтобы задать пользовательское сопоставление, определите функцию со спецификацией операционной точки, opspec, в качестве входного аргумента и выходных аргументов:

Для получения индексов состояния, ввода и вывода на основе путей блоков и имен состояний используйте getStateIndex, getInputIndex, и getOutputIndex. Использование этих команд устойчиво к будущим изменениям модели, таким как добавление состояний модели. Можно также вручную указать индексы. Для получения дополнительной информации о формате indx, indu, и indy, см. getStateIndex, getInputIndex, и getOutputIndex.

Если нет состояний, входов или выходов, используемых пользовательскими функциями ограничения и цели, возвращайте соответствующий выходной аргумент как [].

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

function [indx,indu,indy] = myMapping(opspec)
    indx = [getStateIndex(opspec,'scdTanks/Tank1');
            getStateIndex(opspec,'scdTanks/Tank2');
            getStateIndex(opspec,'scdTanks/Tank3')];
    indu = [];
    indy = [];
end

Добавьте пользовательское сопоставление в спецификацию рабочей точки.

opspec.CustomMappingFcn = @myMapping;

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

function [c_ineq,c_eq] =  myConstraintsMap(x,u,y)
    c_ineq = [];
    c_eq = [x(1)-x(2);     % Tank1 pressure - Tank2 pressure
            x(2)-x(3)+1];  % Tank2 pressure - Tank3 pressure + 1
end

function F =  myObjectiveMap(x,u,y)
    F = max(x(3)-20, 0) + max(16-x(3), 0);
end

Здесь, x, u, и y являются векторами отображенных состояний, входов и выходов соответственно. Эти векторы содержат сопоставленные значения, указанные в indx, indu, и indyсоответственно.

Добавьте обновленные пользовательские функции в спецификацию рабочей точки.

opspec.CustomConstrFcn = @myConstraintsMap;
opspec.CustomObjFcn = @myObjectiveMap;

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

[op3,rpt3] = findop(mdl,opspec,opt);
op3.States

(1.) scdTanks/Inertia
      x: 0            
(2.) scdTanks/Tank1
      x: 15           
(3.) scdTanks/Tank2
      x: 15           
(4.) scdTanks/Tank3
      x: 16           

Добавление аналитических градиентов в пользовательские функции

Для более быстрых или надежных вычислений можно добавить аналитические градиенты к пользовательским ограничениям и объективным функциям. Добавление градиентов может уменьшить число вызовов функций во время оптимизации и потенциально повысить точность результата оптимизации. Если заданы градиенты, их необходимо задать как для пользовательских, так и для целевых функций. (Градиенты для пользовательской обрезки не поддерживаются для Simscape™ моделей.)

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

F = (u(1)+3)^2 + y(1)^2

затем градиент F в отношении u(1) является

G = 2*(u(1)+3)

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

Каждый столбец G_ineq и G_eq содержит градиенты для одного ограничения, а порядок столбцов соответствует порядку строк в соответствующем векторе ограничения. Количество строк в обоих G_ineq и G_eq равно общему количеству состояний, входов и выходов в x, u, и y. Каждый столбец содержит градиенты относительно состояний в x, за которыми следуют входные данные в u, то выходы в y.

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

function [c_ineq,c_eq,G_ineq,G_eq] =  myConstraintsGrad(x,u,y)
    c_ineq = [];
    c_eq = [x(1)-x(2);     % Tank1 pressure - Tank2 pressure
            x(2)-x(3)+1];  % Tank2 pressure - Tank3 pressure + 1
    
    G_ineq = [];
    G_eq = [1  0;
            -1 1;
            0 -1]; 
end

В этой функции строка i из G_eq содержит градиенты относительно состояния x(i).

Аналогично, чтобы добавить градиенты в пользовательскую целевую функцию, укажите дополнительный выходной аргумент G, который содержит градиенты F. G возвращается как вектор столбца с тем же форматом, что и столбцы G_ineq и G_eq.

function [F,G] = myObjectiveGrad(x,u,y)
    F = max(x(3)-20, 0) + max(16-x(3), 0);
    
    if x(3) >= 20
        G = [0 0 1]';
    elseif x(3) <= 16
        G = [0 0 -1]';
    else
        G = [0 0 0]';
    end
end

Поскольку целевая функция в этом примере является кусочно дифференцируемой, значение G зависит от величины давления в Tank3.

Добавьте обновленные пользовательские функции в спецификацию рабочей точки.

opspec.CustomConstrFcn = @myConstraintsGrad;
opspec.CustomObjFcn = @myObjectiveGrad;

Чтобы включить градиенты в алгоритме оптимизации, включите Jacobian вариант оптимизации.

opt.OptimizationOptions.Jacobian = 'on';

Чтобы использовать аналитические якобинцы при обрезке моделей с помощью Диспетчера устойчивого состояния (Steady State Manager) или Линеаризатора модели (Model Linearizer), выберите опцию Включить аналитическую обрезку якобиана (Enable analytic Jacobian trimming).

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

[op4,rpt4] = findop(mdl,opspec,opt);
op4.States

(1.) scdTanks/Inertia
      x: 0            
(2.) scdTanks/Tank1
      x: 15           
(3.) scdTanks/Tank2
      x: 15           
(4.) scdTanks/Tank3
      x: 16           

Результат оптимизации тот же, что и результат для неградирующего решения.

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

rpt3.OptimizationOutput.funcCount

ans =

    25

rpt4.OptimizationOutput.funcCount

ans =

     5

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