Discrete-Valued Variables in Response Optimization (Code)

This example shows how to use response optimization to tune discrete-valued variables. Discrete-valued variables represent model parameters that are restricted to a finite set of values, instead of continuously varying. To use discrete-valued variables for response optimization or parameter estimation, you must use the surrogateopt solver for mixed-integer optimization.

Open the Model

In the sdoMotorPosition model, a PI controller enables the angular position of a DC motor to match a desired reference value. The load on the motor is subject to disturbances, and the controller must reject these disturbances.

open_system('sdoMotorPosition')

Within the Controller subsystem, the PI gains are set by a 12 V source divided down by resistors R1, R2, R3, and R4, as shown in the following diagram.

The proportional gain, Kp, is given by 12*R2/(R1+R2) and the integral gain, Ki, is given by 12*R4/(R3+R4). The initial values of the resistances are R1 = 47 kΩ, R2 = 180 kΩ, R3 = 10 kΩ, and R4 = 10 kΩ.

Specify Discrete Design Variables

Those four resistor values are the model parameters to tune for the optimization. Because resistors are only available in discrete values, use discrete parameters to represent them. To do so, use third input argument to sdo.getParameterFromModel to specify discrete parameters.

DesignVars = sdo.getParameterFromModel('sdoMotorPosition',[],{'R1','R2','R3','R4'});

sdo.getParameterFromModel returns an array of param.Discrete parameter objects that correspond to the model parameters in the order you specified. Set the ValueSet property of each parameter to the set of discrete values allowed for the parameter. Set the Value property to the current value, which must be present in ValueSet.

% R1 values
DesignVars(1).ValueSet = [39 43 47 51 56] * 1e3;
DesignVars(1).Value    = 47*1e3;
% R2 values
DesignVars(2).ValueSet = [150 160 180 200 220] * 1e3;
DesignVars(2).Value    = 180*1e3;
% R3 values
DesignVars(3).ValueSet = [8.2 9.1 10 11 12] * 1e3;
DesignVars(3).Value    = 10*1e3;
% R4 values
DesignVars(4).ValueSet = [8.2 9.1 10 11 12] * 1e3;
DesignVars(4).Value    = 10*1e3;

Specify Design Requirements and Signals

The model applies a step disturbance at 1 second. With the initial resistance values specified in the model, the disturbance causes the motor to deviate by about 20°. The response then settles back to within ±5° of the reference position by 4 seconds after the disturbance. For this example, find new resistor values to improve this specification by 10%, so that the motor deviates no more than 18°, and settles back to within ±4.5° degrees of the reference position by 4 seconds after the disturbance. First, specify these new design requirements using sdo.requirements.SignalBound.

Requirements = struct;
Requirements.UpperBound = sdo.requirements.SignalBound(...
    'BoundMagnitudes', [18 18 ; 4.5 4.5], ...
    'BoundTimes', [1 5 ; 5 15]);
Requirements.LowerBound = sdo.requirements.SignalBound(...
    'BoundMagnitudes', [-4.5 -4.5], ...
    'BoundTimes', [5 15], ...
    'Type', '>=');

To visualize the initial performance against the requirement, first specify signals to log during simulation. In particular, log the output of the Integrator block, which is the angular position of the model, the signal you want to optimize. To do so, create a simulation scenario with sdo.SimulationTest and configure its LoggingInfo property. Also, to prepare the scenario for later use in the optimization, assign the design variables you created as its Parameters property.

Simulator = sdo.SimulationTest('sdoMotorPosition');

Sig_Info = Simulink.SimulationData.SignalLoggingInfo;
Sig_Info.BlockPath = 'sdoMotorPosition/Integrator';
Sig_Info.LoggingInfo.LoggingName = 'Sig';
Sig_Info.LoggingInfo.NameMode = 1;

Simulator.LoggingInfo.Signals = Sig_Info;

Simulator.Parameters = DesignVars;

Now create a plot of the disturbance-rejection requirements. Then, simulate the model and add the logged signal to the plot.

hRequirement = plot( ...
    Requirements.LowerBound.BoundTimes', ...
    Requirements.LowerBound.BoundMagnitudes', ...
    'color','black');
hold on
plot( ...
    Requirements.UpperBound.BoundTimes', ...
    Requirements.UpperBound.BoundMagnitudes', ...
    'color','black');
% Simulator.Parameters = DesignVars;
Simulator = sim(Simulator);
SimLog = find(Simulator.LoggedData, ...
    get_param('sdoMotorPosition','SignalLoggingName'));
Sig_Log = find(SimLog,'Sig');
clrs = colororder;
hOriginal = plot(Sig_Log.Values,'Color',clrs(2,:));
legend([hRequirement hOriginal],'Requirement','Original Variables');

The plot shows that the current resistor values do not quite satisfy the more stringent requirement.

Set up Optimization

The optimization process runs the model many times with different values for the resistor variables. To expedite this, configure the simulator for fast restart.

Simulator = setup(Simulator,FastRestart='on');

Optimization requires an objective or cost function that runs the model and determines whether the disturbance rejection requirements are satisfied. For this example, use the function sdoMotorPosition_optFcn, which is included at the end of the example. This function is called many times by the optimization solver. To pass the function to the optimizer, use an anonymous function with one argument that calls sdoMotorPosition_optFcn.

optimfcn = @(P) sdoMotorPosition_optFcn(P,Simulator,Requirements);

Specify optimization options. Set the optimization method to surrogateopt solver, which is the method that supports tuning of discrete parameters.

Options = sdo.OptimizeOptions;
Options.Method = 'surrogateopt';
Options.MethodOptions.MaxFunctionEvaluations = 100;
Options.MethodOptions.ObjectiveLimit = 0.001;
Options.OptimizedModel = Simulator;

Optimize the Design

To find resistance values optimized for the requirements, call sdo.optimize with the cost function handle, parameters to optimize, and options.

rng('default');   % for reproducibility
[Optimized_DesignVars,Info] = sdo.optimize(optimfcn,DesignVars,Options);

 Optimization started 2024-Sep-05, 19:17:05

                                       Current 
 F-count     max constraint     max constraint     Trial type
       1          0.0759805          0.0759805     initial
       2          0.0751287          0.0751287     random
       3          0.0751287           0.179179     random
       4          0.0447169          0.0447169     random
       5          0.0317818          0.0317818     random
       6          0.0317818           0.261353     random
       7          0.0257711          0.0257711     random
       8          0.0257711          0.0930365     random
       9          0.0257711          0.0971938     random
      10          0.0257711           0.200359     random
      11          0.0062997          0.0062997     random
      12          0.0062997           0.206374     random
      13          0.0062997          0.0818387     random
      14          0.0062997          0.0585416     random
      15          0.0062997           0.114002     random
      16          0.0062997           0.112485     random
      17          0.0062997          0.0759805     random
      18          0.0062997          0.0709489     random
      19          0.0062997           0.280136     random
      20          0.0062997           0.099985     random
      21        -0.00759693        -0.00759693     adaptive
      22        -0.00759693          0.0691728     random
      23        -0.00759693          0.0347164     random
      24        -0.00759693          0.0977682     random
      25        -0.00759693          0.0589627     random
      26        -0.00759693           0.131778     random
      27        -0.00759693           0.147347     random
      28        -0.00759693           0.104488     random
      29        -0.00759693          0.0534262     random
      30        -0.00759693             0.1643     random
      31        -0.00759693          0.0419731     random
      32        -0.00759693          0.0396579     random
      33        -0.00759693           0.077945     random
      34        -0.00759693          0.0883493     random
      35        -0.00762259        -0.00762259     random
      36        -0.00762259          0.0455853     random
      37        -0.00762259            0.13976     random
      38        -0.00762259           0.140634     random
      39        -0.00762259          0.0662368     random
      40        -0.00762259           0.256738     random
      41        -0.00762259          0.0348764     random
      42        -0.00762259            0.18995     random
      43        -0.00762259          0.0615952     random
      44        -0.00762259          0.0597358     random
      45        -0.00762259           0.015983     random
      46        -0.00762259          0.0770592     random
      47        -0.00762259           0.174666     random
      48        -0.00762259           0.081758     random
      49        -0.00762259          0.0240828     random
      50        -0.00762259           0.173292     random

                                       Current 
 F-count     max constraint     max constraint     Trial type
      51        -0.00762259           0.116519     random
      52        -0.00762259           0.217099     random
      53        -0.00762259           0.104371     random
      54        -0.00762259          0.0937562     random
      55        -0.00762259           0.104225     random
      56        -0.00762259          0.0380129     random
      57        -0.00762259           0.150819     random
      58        -0.00762259          0.0636666     random
      59        -0.00762259          0.0320499     random
      60        -0.00762259          0.0273158     random
      61        -0.00762259          0.0798937     random
      62        -0.00762259           0.174146     random
      63        -0.00762259          0.0166127     random
      64        -0.00762259         0.00433873     adaptive
      65        -0.00762259           0.294106     random
      66        -0.00518218        -0.00518218     adaptive
      67        -0.00518218          0.0319864     random
      68        -0.00518218           0.266843     random
      69        -0.00518218          0.0342618     random
      70        -0.00518218          0.0982987     random
      71        -0.00518218          0.0419731     random
      72        -0.00518218          0.0396579     random
      73        -0.00518218           0.299815     random
      74        -0.00518218          0.0712971     random
      75        -0.00518218          0.0709489     random
      76        -0.00518218          0.0589627     random
      77        -0.00518218           0.131778     random
      78        -0.00518218         0.00653156     random
      79        -0.00518218           0.211877     random
      80        -0.00518218          0.0967155     random
      81        -0.00518218          0.0689487     random
      82        -0.00518218          0.0257711     random
      83        -0.00518218          0.0930365     random
      84        -0.00518218           0.114002     random
      85        -0.00518218           0.112485     random
      86        -0.00518218            0.06706     random
      87        -0.00518218          0.0751287     random
      88        -0.00737024        -0.00737024     adaptive
      89        -0.00737024          0.0759805     random
      90        -0.00737024            0.10447     random
      91        -0.00737024           0.190644     random
      92        -0.00737024           0.142755     random
      93        -0.00737024           0.104193     random
      94        -0.00737024           0.100237     random
      95        -0.00737024           0.097787     random
      96        -0.00737024          0.0735855     random
      97        -0.00737024           0.137557     random
      98        -0.00737024          0.0783709     random
      99        -0.00737024          0.0875674     random
     100        -0.00737024          0.0794652     random
     101        -0.00737024        -0.00737024     best value
The current solution is feasible. surrogateopt stopped because it exceeded the function evaluation limit set by the 'MethodOptions.MaxFunctionEvaluations' property in the sdo.OptimizeOptions object.
If the solution needs to be improved, you could try increasing the function evaluation limit.

Evaluate Optimized Design

Plot the model response after optimization.

Simulator.Parameters = Optimized_DesignVars;
Simulator = sim(Simulator);
SimLog = find(Simulator.LoggedData, ...
    get_param('sdoMotorPosition','SignalLoggingName'));
Sig_Log = find(SimLog,'Sig');
clrs = colororder;
hOptimized = plot(Sig_Log.Values, 'Color', clrs(1,:));
legend([hRequirement hOriginal hOptimized], ...
    'Requirement','Original Variables','Optimized Variables');

Now the motor position is within the spec of the disturbance rejection requirement.

Update Model with Optimized Parameter Values

sdo.Optimize returns tuned versions of the parameters in the array Optimized_DesignVars. Each entry in this array is a param.Discrete parameter object whose Value property is set to the optimized value. For instance, examine the new value of R2.

Optimized_DesignVars(2).Value

ans = 
220000

Update the model with the optimized parameter values.

sdo.setValueInModel('sdoMotorPosition',Optimized_DesignVars);

Restore the simulator fast restart settings.

Simulator = restore(Simulator);

Conclusion

In this example you tuned variables to improve the disturbance rejection characteristics of a motor controller. The variables being tuned were electrical component parts that could only take on discrete values, rather than continuous values. You used sdo.getParameterFromModel to specify the variables as discrete, and used the surrogateopt solver to tune the parameters.

Objective Function

The function sdoMotorPosition_optFcn is called at each iteration of the optimization problem with a set of parameter values, P. The function returns the objective value and constraint violations, Vals, to the optimization solver. See sdo.optimize for a more detailed description of the function signature.

function Vals = sdoMotorPosition_optFcn(P,Simulator,Requirements)
%SDOMOTORPOSITION_OPTFCN
%

% Simulate the model.
Simulator.Parameters = P;
Simulator = sim(Simulator);

% Retrieve logged signal data.
SimLog = find(Simulator.LoggedData, ...
    get_param('sdoMotorPosition','SignalLoggingName'));
Sig_Log = find(SimLog,'Sig');

% Evaluate the design requirements.
Cleq_UpperBound = evalRequirement( ...
    Requirements.UpperBound,Sig_Log.Values);
Cleq_LowerBound = evalRequirement( ...
    Requirements.LowerBound,Sig_Log.Values);

% Collect the evaluated design requirement values in a structure to
% return to the optimization solver.
Vals.Cleq = [...
    Cleq_UpperBound(:); ...
    Cleq_LowerBound(:)];
end