Модель 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

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