Estimate Model Parameters of a Symbolically Derived Plant Model in Simulink

This example uses Simulink Design Optimization™ to estimate the unknown capacitance and initial voltage of a symbolically derived algebraic model of a simple resistor-capacitor (RC) circuit. The example solves the same problem and uses the same experimental data as Estimate Model Parameters and Initial States (Code) (Simulink Design Optimization), but uses the closed-form solution for the RC circuit instead of the differential form.

This example uses Symbolic Math Toolbox™ capabilities to:

• Solve ordinary differential equations (ODE) using dsolve

• Convert an analytical result into a Simulink® block using matlabFunctionBlock

You perform design optimization to estimate the capacitance and initial voltage values of the analytical RC circuit. In particular, you match the experimental output voltage values with the simulated values.

Solve Equation for RC Circuit

Define and solve the following differential equation for the RC circuit.

Here, v2(t) is the output voltage across capacitor C1, v1 is the constant voltage across resistor R1, and v20 is the initial voltage across the capacitor. Use dsolve to solve the equation.

syms C1 R1 v1 v20 real 
syms v2(t)
deq = (v1 - v2)/R1 - C1*diff(v2,t);
v2sol = dsolve(deq, v2(0) == v20)

v2sol = 
$$v_1-{\mathrm{e}}^{-\frac{t}{C_1\hspace{0.17em}{R}_1}}\hspace{0.17em}\left(v_1-v_{20}\right)$$

Use subs to evaluate the solution for an R1 value of 10 kOhm and v1 value of 5 V.

v2sol = vpa(subs(v2sol,[R1,v1],[10e3,5]))

v2sol = 
$${\mathrm{e}}^{-\frac{0.0001\hspace{0.17em}t}{C_1}}\hspace{0.17em}\left(v_{20}-5.0\right)+5.0$$

Create Model with Block Representing RC Circuit

First, create a new Simulink model.

myModel = 'rcSymbolic';
new_system(myModel);
load_system(myModel);

Use matlabFunctionBlock to convert the symbolic result for the output voltage to a Simulink block representing the RC plant model. matlabFunctionBlock adds this new block to the model.

blockName = 'closedFormRC_block';
rcBlock = strcat(myModel,'/',blockName);
myVars = [C1,v20,t];
matlabFunctionBlock(rcBlock,v2sol,...
    'vars',myVars,...
    'functionName','myRC',...
    'outputs',{'v2'});

Add More Blocks

Add and arrange other blocks with positions and sizes relative to the RC block.

rcBlockPosition = get_param(rcBlock,'position');
rcBlockWidth = rcBlockPosition(3)-rcBlockPosition(1);
rcBlockHeight = rcBlockPosition(4)-rcBlockPosition(2);
constantBlock = 'built-in/Constant';
timeBlock = 'simulink/Sources/Ramp';
outputBlock = 'built-in/Outport';

C1 and v20 are the parameters to estimate. First, introduce and initialize these parameters in the MATLAB® workspace, using the initial values of 460 $$\mu F$$ and 1 V, respectively. Then create constant blocks for both parameters.

C1val = 460e-6;
v20val = 1.0;
posX = rcBlockPosition(1)-rcBlockWidth*2;
posY = rcBlockPosition(2)-rcBlockHeight*3/4;
pos = [posX,posY,posX+rcBlockWidth/2,posY+rcBlockHeight/2];
add_block(constantBlock,strcat(myModel,'/C1'),'Value','C1val',...
    'Position',pos);
pos = pos + [0 rcBlockHeight 0 rcBlockHeight];
add_block(constantBlock,strcat(myModel,'/v20'),'Value','v20val',...
    'Position',pos);

Add a ramp for time.

pos = pos + [0 rcBlockHeight 0 rcBlockHeight];
add_block(timeBlock,strcat(myModel,'/t'),'Slope','1','Position',pos);

Add an output port.

pos = [rcBlockPosition(1)+2*rcBlockWidth,...
    rcBlockPosition(2)+rcBlockHeight/4,...
    rcBlockPosition(1)+2*rcBlockWidth+rcBlockWidth/2,...
    rcBlockPosition(2)+rcBlockHeight/4+rcBlockHeight/2];
add_block(outputBlock,strcat(myModel,'/v2'),'Port','1','Position',pos);

Now, wire blocks in the model. The model is ready for Simulink Design Optimization.

myAddLine = @(k) add_line(myModel,...
    strcat(char(myVars(k)),'/1'),...
    strcat(blockName,'/',num2str(k)),...
    'autorouting','on');
arrayfun(myAddLine,(1:numel(myVars)));
add_line(myModel,strcat(blockName,'/1'),'v2/1','autorouting','on');
open_system(myModel);

Estimate Parameters

Get the measured data.

load sdoRCCircuit_DataFromExperiment

The variables time and data are loaded into the workspace. data is the measured capacitor voltage for times time.

Create an sdo.Experiment object to store the experimental voltage data.

Exp = sdo.Experiment(myModel);

Create an object to store the measured capacitor voltage output.

Voltage = Simulink.SimulationData.Signal;
Voltage.Name      = 'Voltage';
Voltage.BlockPath = rcBlock;
Voltage.PortType  = 'outport';
Voltage.PortIndex = 1;
Voltage.Values    = timeseries(data,time);

Add the measured capacitor data to the experiment as the expected output data.

Exp.OutputData = Voltage;

Get the parameters. Set a minimum value for C1. Note that you already specified the initial guesses.

c1param = sdo.getParameterFromModel(myModel,'C1val');
c1param.Minimum = 0;
v20param = sdo.getParameterFromModel(myModel,'v20val');

Define the objective function for estimation. The code for sdoRCSymbolic_Objective used in this example is given in the Helper Functions section at the end of the example.

estFcn = @(v) sdoRCSymbolic_Objective(v,Exp,myModel);

Collect the model parameters to be estimated.

v = [c1param;v20param];

Because the model is entirely algebraic, turn off warning messages that instruct you to select a discrete solver.

set_param(myModel,'SolverPrmCheckMsg','none');

Estimate the parameters.

opt = sdo.OptimizeOptions;
opt.Method = 'lsqnonlin';
vOpt = sdo.optimize(estFcn,v,opt);

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

                                          First-order 
 Iter F-count        f(x)      Step-size  optimality
    0      5      27.7093            1
    1     10      2.86889        1.919         2.94
    2     15      1.53851       0.3832        0.523
    3     20      1.35137       0.3347        0.505
    4     25      1.34473      0.01374      0.00842
    5     30      1.34472     0.002686      0.00141
Local minimum possible.

lsqnonlin stopped because the final change in the sum of squares relative to 
its initial value is less than the value of the function tolerance.

Show the estimated values.

fprintf('C1 = %e v20 = %e\n',vOpt(1).Value, vOpt(2).Value);

C1 = 2.261442e-04 v20 = 2.359446e+00

Compare Simulated and Experimental Data

Update the experiment with the estimated capacitance and capacitor initial voltage values.

Exp = setEstimatedValues(Exp,vOpt);

Simulate the model with the estimated parameter values and compare the simulated output with the experimental data.

Simulator = createSimulator(Exp);
Simulator = sim(Simulator);
SimLog    = find(Simulator.LoggedData,get_param(myModel,'SignalLoggingName'));
Voltage   = find(SimLog,'Voltage');
plot(time,data,'ro',Voltage.Values.Time,Voltage.Values.Data,'b')
title('Simulated and Measured Responses')
legend('Measured Voltage','Simulated Voltage','Location','Best')

close_system(myModel,0);

Helper Functions

function vals = sdoRCSymbolic_Objective(v,Exp,myModel) 
    r = sdo.requirements.SignalTracking;
    r.Type      = '==';
    r.Method    = 'Residuals';
    r.Normalize = 'off';
    Exp  = setEstimatedValues(Exp,v);
    Simulator = createSimulator(Exp);
    Simulator = sim(Simulator);
    SimLog  = find(Simulator.LoggedData,get_param(myModel,'SignalLoggingName'));
    Voltage = find(SimLog,'Voltage');
    VoltageError = evalRequirement(r,Voltage.Values,Exp.OutputData(1).Values);
    vals.F = VoltageError(:);
end