# Parameter Estimation with Constraints

Since R2023b

This example shows how to perform parameter estimation while also imposing constraints the model needs to obey. In this example, you estimate the parameters of an engine throttle system.

### Engine Throttle Model

The throttle controls the air flow into the intake manifold of an engine. The throttle body contains a valve that opens when the driver presses down on the accelerator pedal. This lets more air enter the cylinders.

A DC motor controls the opening angle of the butterfly valve. There is also a spring attached to the valve to return it to its closed position when the DC motor is de-energized. The amount of rotation of the valve is limited to approximately 90 degrees. Therefore, if a large command input is applied to the motor, the valve hits the hard stops preventing it from rotating further.

The DC motor is modeled as a torque gain and a time-delay input with parameters `Kt` and `input_delay`. The butterfly valve is modeled as a mass-spring-damper system with parameters `J`, `c` and `k`. This system is augmented with hard stops to limit the valve opening to 90 degrees. The model components are known, however the parameter values of the system are not known accurately and need to be estimated. Additionally, the torque produced by the motor is bounded. Therefore, you need to impose this constraint on the parameter estimation process.

Open the model.

`open_system('spe_engine_throttle')`

### Specify Model Parameters to Estimate

The parameters to estimate are:

• Motor delay `input_delay`.

• Throttle moment of inertia `J`, damping `c`, and spring constant, `k`.

`p = sdo.getParameterFromModel('spe_engine_throttle',{'J','c','input_delay','k'});`

Set initial guesses for these parameters.

```p(1).Value = 0.2; p(2).Value = 9; p(3).Value = 0.01; p(4).Value = 18;```

Since these are all physical parameters, specify the minimum values as 0 or positive.

```p(1).Minimum = 0; p(2).Minimum = 0; p(3).Minimum = 1e-6; p(3).Maximum = 0.1; p(4).Minimum = 0;```

### Specify Constraints

The motor torque is bounded and this constraint must be satisfied during parameter tuning. Describe this constraint so you can impose it as a requirement during parameter estimation.

```Requirements = struct; boundMagnitude = 15; Requirements.SignalBound = sdo.requirements.SignalBound(... 'BoundTimes', [0 0.5], ... 'BoundMagnitudes', boundMagnitude*[1 1], ... 'Type', '<='); Requirements.SignalBound2 = sdo.requirements.SignalBound(... 'BoundTimes', [0 0.5], ... 'BoundMagnitudes', -boundMagnitude*[1 1], ... 'Type', '>=');```

To evaluate the requirement, log the model signals on which these requirements are based. Specify the model signals to log during model simulation.

```Sig_Info = Simulink.SimulationData.SignalLoggingInfo; Sig_Info.BlockPath = 'spe_engine_throttle/Throttle/Velocity'; Sig_Info.LoggingInfo.LoggingName = 'Sig'; Sig_Info.LoggingInfo.NameMode = 1;```

Define a simulator, which is runs the model with different parameter values. Add the signal logging specification to the simulator.

```Simulator = sdo.SimulationTest('spe_engine_throttle'); Simulator.LoggingInfo.Signals = Sig_Info;```

### Define Estimation Experiments

The estimation experiment compares the measured data with the model output, at the corresponding model signal. Define an experiment to associate measured data with the model signal.

```Time = time1; MotorCommand = input1; ThrottlePosition = position1; EstimationData = sdo.Experiment('spe_engine_throttle');```

Specify the measured input data and associate it with the corresponding model input.

```EstimationData_Sig_Input = Simulink.SimulationData.Signal; EstimationData_Sig_Input.Values = timeseries(MotorCommand,Time); EstimationData_Sig_Input.BlockPath = 'spe_engine_throttle/Input'; EstimationData_Sig_Input.PortType = 'inport'; EstimationData_Sig_Input.PortIndex = 1; EstimationData_Sig_Input.Name = 'spe_engine_throttle/Input:1'; EstimationData.InputData = EstimationData_Sig_Input;```

Specify the measured output data and associate it with the corresponding model signal.

```EstimationData_Sig_Output = Simulink.SimulationData.Signal; EstimationData_Sig_Output.Values = timeseries(ThrottlePosition,Time); EstimationData_Sig_Output.BlockPath = 'spe_engine_throttle/Throttle'; EstimationData_Sig_Output.PortType = 'outport'; EstimationData_Sig_Output.PortIndex = 1; EstimationData_Sig_Output.Name = 'spe_engine_throttle/Throttle:1'; EstimationData.OutputData = EstimationData_Sig_Output;```

Add the experiment to the model simulator, so the model simulation uses the information in the experiment.

`Simulator = createSimulator(EstimationData,Simulator);`

Try the nominal parameter values in the model. To run the model and plot results, use the function `spe_engine_throttle_Plots`.

`spe_engine_throttle_Plots(p,Simulator,EstimationData,Requirements);`

In these plots, the fit between measured data and simulated throttle position is fairly good, but the motor torque does not satisfy the requirement bound, indicating that you can estimate better parameter values.

### Optimization Options

To obtain better parameter values, set up parameter estimation. Specify optimization options and use the `lsqnonlin` optimization solver.

```Options = sdo.OptimizeOptions; Options.Method = 'lsqnonlin'; Options.MethodOptions.ConstraintTolerance = 1e-06; Options.MethodOptions.FunctionTolerance = 1e-06; Options.MethodOptions.OptimalityTolerance = 1e-06; Options.MethodOptions.StepTolerance = 1e-06; Options.OptimizedModel = Simulator;```

### Create Estimation Objective Function

Create a function that the optimizer calls at each optimization iteration to compute the estimation cost. Use an anonymous function with one argument that calls `spe_engine_throttle_optFcn`.

`optimfcn = @(P) spe_engine_throttle_optFcn(P,Simulator,EstimationData,Requirements);`

### Estimate the Parameters

Call `sdo.optimize` with the estimation objective function handle, parameters to estimate, and options.

`[pOpt,Info] = sdo.optimize(optimfcn,p,Options);`
``` Optimization started 2024-Feb-13, 01:27:18 First-order Iter F-count f(x) constraint Step-size optimality 0 13 0.0574662 0.129775 0 1 27 0.0543575 0.082508 0.04755 11.1 2 41 0.100572 0.018302 0.1103 9.41 3 54 0.139486 0 0.06346 3.08 4 67 0.166934 0 0.02637 0.1 5 80 0.106579 0 0.08005 1.47 6 93 0.0751622 0 0.06821 0.332 7 106 0.0697339 0.000121072 0.0336 0.356 8 120 0.067084 0 0.07122 0.139 9 133 0.0656637 0 0.05265 0.0764 10 146 0.0650715 0 0.03425 0.137 11 159 0.0648462 0 0.02148 0.188 12 173 0.0647772 0 0.01659 0.313 13 186 0.0647238 0 0.005592 0.231 14 203 0.064699 0 0.000388 0.234 15 230 0.0646956 0 2.707e-06 0.194 Local minimum possible. Constraints satisfied. lsqnonlin stopped because the size of the current step is less than the value of the step size tolerance and constraints are satisfied to within the value of the constraint tolerance. ```

The progress display provides information about the iterations from parameter estimation. By the end of the estimation process, the fit between measured data and simulated throttle position is fairly good, as indicated by a small value in the `f(x)` column. In addition, the requirement is satisfied, as indicated by 0 or a negative value in the `constraint` column.

Try the optimized parameter values in the model.

`spe_engine_throttle_Plots(pOpt,Simulator,EstimationData,Requirements);`

The plots show a good fit between measured data and simulated throttle position, and the motor torque satisfies the requirement bound.

Update the experiment with the estimated parameter values.

`EstimationData = setEstimatedValues(EstimationData,pOpt);`

Update the model with the optimized parameter values.

`sdo.setValueInModel('spe_engine_throttle',pOpt);`

In summary, the parameters of the engine throttle model are estimated to make the model output match measured data, while also imposing a constraint on the simulated motor torque, to keep it consistent with a known bound.

Close the model.

`close_system("spe_engine_throttle");`

### Example Helper Functions

`type spe_engine_throttle_Plots`
```function spe_engine_throttle_Plots(P,Simulator,EstimationData,Requirements) %% Evaluate Requirements % Update parameter values, and simulate the model Simulator.Parameters = P; Simulator = sim(Simulator); % Retrieve logged signal data SimLog = find(Simulator.LoggedData,get_param('spe_engine_throttle','SignalLoggingName')); Sig_Log = find(SimLog,'Sig'); %% Evaluate Experiments % Update parameter values in the experiment, and simulate the model EstimationData = setEstimatedValues(EstimationData,P); Simulator = createSimulator(EstimationData,Simulator); strOT = mat2str(EstimationData.OutputData(1).Values.Time); Simulator = sim(Simulator, 'OutputOption', 'AdditionalOutputTimes', 'OutputTimes', strOT); % Retrieve logged signal data SimLog = find(Simulator.LoggedData,get_param('spe_engine_throttle','SignalLoggingName')); Sig = find(SimLog,EstimationData.OutputData.Name); %% Plot Fit to Measured Data in Experiment colorOrder = colororder('default'); subplot(1,2,1); plot(EstimationData.OutputData.Values.Time, EstimationData.OutputData.Values.Data, ... 'Color', colorOrder(1,:), 'DisplayName', 'Measured'); line(Sig.Values.Time, Sig.Values.Data, ... 'Color', colorOrder(2,:), 'DisplayName', 'Simulation'); legend('Location','southeast'); title('Experiment'); ylabel('Throttle Position'); xlabel('Time'); %% Plot Requirement subplot(1,2,2); plot(Sig_Log.Values.Time, Sig_Log.Values.Data, 'Color', colorOrder(5,:)); line(Requirements.SignalBound.BoundTimes, Requirements.SignalBound.BoundMagnitudes, ... 'Color', colorOrder(3,:), 'LineWidth', 1.25); if numel(fieldnames(Requirements)) > 1 line(Requirements.SignalBound2.BoundTimes, Requirements.SignalBound2.BoundMagnitudes, ... 'Color', colorOrder(3,:), 'LineWidth', 1.25); end legend('Model Signal', 'Bound', 'Location','southeast'); title('Requirement'); ylabel('Motor Torque'); xlabel('Time'); end ```
`type spe_engine_throttle_optFcn`
```function Vals = spe_engine_throttle_optFcn(P,Simulator,EstimationData,Requirements) %SPE_ENGINE_THROTTLE_OPTFCN % % Function called at each iteration of the estimation problem. % % The function is called with a set of parameter values, P, and returns % the estimation cost, Vals, to the optimization solver. % % See the sdoExampleCostFunction function and sdo.optimize for a more % detailed description of the function signature. % %% Evaluate Requirements % % Simulate the model. Simulator.Parameters = P; Simulator = sim(Simulator); % Retrieve logged signal data. SimLog = find(Simulator.LoggedData,get_param('spe_engine_throttle','SignalLoggingName')); Sig_Log = find(SimLog,'Sig'); % Evaluate the design requirements. Cleq_SignalBound = evalRequirement(Requirements.SignalBound, Sig_Log.Values); Cleq_SignalBound2 = evalRequirement(Requirements.SignalBound2, Sig_Log.Values); %% Evaluate Experiments % % Define a signal tracking requirement to compute how well the model % output matches the experiment data. r = sdo.requirements.SignalTracking('Method','Residuals'); %% % Update the experiment with the estimated parameter values. EstimationData = setEstimatedValues(EstimationData,P); %% % Simulate the model and compare model outputs with measured experiment % data. F_r = []; Simulator = createSimulator(EstimationData,Simulator); strOT = mat2str(EstimationData.OutputData(1).Values.Time); Simulator = sim(Simulator, 'OutputOption', 'AdditionalOutputTimes', 'OutputTimes', strOT); SimLog = find(Simulator.LoggedData,get_param('spe_engine_throttle','SignalLoggingName')); Sig = find(SimLog,EstimationData.OutputData.Name); Error = evalRequirement(r,Sig.Values,EstimationData.OutputData.Values); F_r = [F_r; Error(:)]; %% Return Values. % % Return the evaluated estimation cost and requirement values in a % structure to the optimization solver. %Return values from cost function %to optimization solver Vals.F = F_r; Vals.Cleq = [... Cleq_SignalBound(:) ; ... Cleq_SignalBound2(:) ... ]; end ```