## Compute Operating Points Using Custom Constraints and Objective Functions

Typically, when computing a steady-state operating point for a Simulink® model using an optimization-based search, you specify known fixed values or bounds to constrain your model states, inputs, or outputs. However, some systems or applications require additional flexibility in defining the optimization search parameters.

For such systems, you can specify custom constraints, an additional optimization objective function, or both. When the software computes a steady-state operating point, it applies these custom constraints and objective function in addition to the standard state, input, and output specifications.

You can specify custom equality and inequality constraints as algebraic combinations of model states, inputs, and outputs. These constraints let you limit the operating point search space by specifying known relationships between inputs, outputs, and states. For example, you can specify that one model state is the sum of two other states.

You can also specify a custom scalar objective function as an algebraic combination of model states, inputs, and outputs. Using the objective function you can optimize the steady-state operating point based on your application requirements. For example, suppose that your model has multiple potential equilibrium points. You can specify an objective function to find the steady-state point with the minimum input energy.

For complex models, you can specify a custom mapping function that selects a subset of the model inputs, outputs, and states to pass to the custom cost and constraint functions.

You can specify custom optimization functions when trimming your model:

• At the command line: Create an operating point specification using `operspec`, and specify the custom functions using the `CustomConstrFcn`, `CustomCostFcn`, and `CustomMappingFcn` properties of the specification.

• Using Steady State Solver: On the Specification tab, click Trim Options. In the Trim Options dialog box, in the Custom Optimization Functions section, specify the function names.

• Using the Linear Analysis Tool: On the Linear Analysis tab, in the Operating Point drop-down list, click Trim Model. In the Trim the model dialog box, on the Options tab, in the Custom Optimization Functions section, specify the function names. The following example shows how to create custom optimization functions and how to trim a model at the command line using these custom functions.

For this example, use a model of three tanks connected with each other by orifices.

```mdl = 'scdTanks'; open_system(mdl) ``` The flow between Tank1 and Tank2 is desired. The flow between Tank2 and Tank3 is undesired unavoidable leakage.

At the expected steady state of this system:

• Tank1 and Tank2 have the same pressure.

• Tank2 and Tank3 have an almost constant pressure difference of `1` that compensates a load.

Due to the weak connectivity between Tank1 and Tank2, it is difficult to trim the model such that the pressures in Tank1 and Tank2 are equal.

### Trim Model Without Customizations

Create a default operating point specification for the model. The specification configures all three tank pressures as free states that must be at steady state in the trimmed operating point.

```opspec = operspec(mdl); ```

Create an option set for trimming the model, suppressing the Command Window display of the operating point search report. The specific trimming options depend on your application. For this example, use nonlinear least squares optimization.

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

Trim the model, and view the trimmed tank pressures.

```[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 ```

The trimmed pressures in Tank1 and Tank2 do not match. Thus, the default operating point specification fails to find an operating point that meets the expected steady-state requirements. If you reduce the constraint tolerance, `opt.OptimizationOptions.TolCon`, you cannot achieve a feasible steady-state solution due to the leakage between Tank2 and Tank3.

To specify custom constraints, define a function in the current working folder or on the MATLAB path with input arguments:

• `x` - Operating point specification states, specified as a vector.

• `u` - Operating point specification inputs, specified as a vector.

• `y` - Operating point specification outputs, specified as a vector.

and output arguments:

• `c_ineq` - Inequality constraints which must satisfy `c_ineq <= 0` during trimming, returned as a vector.

• `c_eq` - Equality constraints which must satisfy `c_eq = 0` during trimming, returned as a vector.

Each element of `c_ineq` and `c_eq` specifies a single constraint. Define the specific constraints for your application as algebraic combinations of the states, inputs, and outputs. If there are no custom equality or inequality constraints, return the corresponding output argument as `[]`.

For this example, to satisfy the conditions of the expected steady state, define the following custom constraint function.

```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 ```

The first entry of `c_eq` constrains the pressures of Tank1 and Tank2 to be the same value. The second equality constraint defines the pressure drop between Tank2 and Tank3.

Add the custom constraint function to the operating point specification.

```opspec.CustomConstrFcn = @myConstraints; ```

Trim the model using the revised operating point specification that contains the custom constraints, and view the trimmed state values.

```[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 ```

Trimming the model with the custom constraint function produces an operating point with equal pressures in the first and second tanks, as expected. Also, as expected, there is a pressure differential of `1` between the third and second tanks.

To examine the final values of the specified constraints, you can check the `CustomEqualityConstr` and `CustomInequalityConstr` properties of the operating point search report.

```rpt1.CustomEqualityConstr ```
```ans = 1.0e-06 * -0.0001 -0.1540 ```

The near-zero values indicate that the equality constraints are satisfied.

To specify a custom objective function, define a function with the same input arguments as the custom constraint function (`x`, `u`, and `y`), and output argument `F`. `F` is an objective function value to be minimized during trimming, returned as a scalar.

Define the objective function for your application as an algebraic combination of the states, inputs, and outputs.

For this example, assume that you want to keep the pressure in Tank3 in the range [16,20]. However, this condition is not always feasible. Therefore, rather than impose hard constraints, add an objective function to incur a penalty if the pressures are not in the [16,20] range. To do so, define the following custom objective function.

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

Add the custom objective function to the operating point specification object.

```opspec.CustomObjFcn = @myObjective; ```

Trim the operating point using both the custom constraints and the custom objective function, and view the trimmed state values.

```[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 ```

In the trimmed operating point, the pressure in Tank3 is within the [16,20] range specified in the custom objective function.

To view the final value of the scalar objective function, check the `CustomObj` property of the operating point search report.

```rpt2.CustomObj ```
```ans = 0 ```

For complex models, you can define a custom mapping that selects a subset of the model states, inputs, and outputs to pass to the custom constraint and objective functions. Doing so simplifies the constraint and objective functions by eliminating unneeded states, inputs, and outputs.

To specify a custom mapping, define a function with your operating point specification, `opspec`, as an input argument, and output arguments:

• `indx` - Indices of mapped states

• `indu` - Indices of mapped inputs

• `indy` - Indices of mapped outputs

To obtain state, input, and output indices based on block paths and state names use `getStateIndex`, `getInputIndex`, and `getOutputIndex`. Using these commands is robust to future model changes, such as the addition of model states. Alternatively, you can manually specify the indices. For more information on the format of `indx`, `indu`, and `indy`, see `getStateIndex`, `getInputIndex`, and `getOutputIndex`.

If there are no states, inputs, or outputs used by the custom constraint and objective functions, return the corresponding output argument as `[]`.

For this example, create a mapping that includes only the pressure states for the three tanks. To do so, define the following custom mapping function.

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

Add the custom mapping to the operating point specification.

```opspec.CustomMappingFcn = @myMapping; ```

When you use a custom mapping function, the indices for the states, inputs, and outputs in your custom constraint and objective functions must be relative to the order specified in the mapping function. Update the custom constraint and objective functions with the new mapping.

```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 ```

Here, `x`, `u`, and `y` are vectors of mapped states, inputs, and outputs, respectively. These vectors contain the mapped values specified in `indx`, `indu`, and `indy`, respectively.

Add the updated custom functions to the operating point specification.

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

Trim the model using the custom mapping, and view the trimmed states, which match the previous results in `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 ```

For faster or more reliable computations, you can add analytic gradients to your custom constraint and objective functions. Adding gradients can reduce the number of function calls during optimization and potentially improve the accuracy of the optimization result. If you specify gradients, you must specify them for both the custom constraint and objective functions. (Gradients for custom trimming are not supported for Simscape™ models.)

To define the gradient of a given constraint or objective function, take the derivative of the function with respect to a given state, input, or output. For example, if the objective function is

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

then the gradient of `F` with respect to `u(1)` is

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

• `G_ineq` - Gradient array for the inequality constraints

• `G_eq` - Gradient array for the equality constraints

Each column of `G_ineq` and `G_eq` contains the gradients for one constraint, and the order of the columns matches the order of the rows in the corresponding constraint vector. The number of rows in both `G_ineq` and `G_eq` is equal to the total number of states, inputs, and outputs in `x`, `u`, and `y`. Each column contains gradients with respect to the states in `x`, followed by the inputs in `u`, then the outputs in `y`.

For this example, add gradients to the constraint function that uses the custom mapping. You do not have to specify a custom mapping when using gradients. However, defining gradients is simpler when using mapped subsets of states, inputs, and outputs.

```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 ```

In this function, row `i` of `G_eq` contains gradients with respect to state `x(i)`.

Similarly, to add gradients to your custom objective function, specify an additional output argument `G`, which contains the gradients of `F`. `G` is returned as a column vector with the same format as the columns of `G_ineq` and `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 ```

Because the objective function in this example is piecewise differentiable, the value of `G` depends on the value of the pressure in Tank3.

Add the updated custom functions to the operating point specification.

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

To enable gradients in the optimization algorithm, enable the `Jacobian` optimization option.

```opt.OptimizationOptions.Jacobian = 'on'; ```

To use analytic Jacobians when trimming models using Steady State Solver or the Linear Analysis Tool, select the Enable analytic Jacobian trimming option. Trim the model using the custom functions with gradients, and view the trimmed states.

```[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 ```

The optimization result is the same as the result for the nongradient solution.

To see if the gradients improve the optimization efficiency, view the operating point search reports. For example, compare the number function evaluations for the solution:

```rpt3.OptimizationOutput.funcCount ```
```ans = 25 ```

```rpt4.OptimizationOutput.funcCount ```
```ans = 5 ```

For this example, adding the analytical gradients decreases the number of function calls during optimization. 