Оцените параметры модели символически выведенной модели объекта управления в Simulink

Этот пример использует Simulink Design Optimization™, чтобы оценить неизвестные емкость и начальное напряжение символьно выведенной алгебраической модели схемы простой резистор-конденсатор (RC). Пример решает ту же задачу и использует те же экспериментальные данные, что и Estimate Параметры модели и Initial States, но использует решение закрытой формы для схемы RC вместо дифференциальной формы.

Этот пример использует возможности Symbolic Math Toolbox™ для:

  • Решить обыкновенные дифференциальные уравнения (ОДУ) используя dsolve

  • Преобразуйте аналитический результат в блок Simulink с помощью matlabFunctionBlock

Вы выполняете проектную оптимизацию, чтобы оценить емкость и начальные значения напряжения аналитической схемы RC. В частности, вы сопоставляете экспериментальные выходы напряжения с моделируемыми значениями.

Решение уравнения для RC-схемы

Определите и решите следующее дифференциальное уравнение для RC-схемы.

Здесь, v2(t) - выход напряжение на конденсаторе C1, v1 - постоянное напряжение на резисторе R1, и v20 - начальное напряжение на конденсаторе. Использование dsolve для решения уравнения.

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

v1-e-tC1R1v1-v20v1 - exp ((-t/( C1 * R1))) * (v1 - v20)

Использование subs для оценки решения для R1 значение 10 кОм и v1 значение 5 В.

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

e-0.0001tC1v20-5.0+5.0exp ((- (vpa ('0.0001') * t )/C1)) * (v20 - vpa ('5.0')) + vpa ('5.0')

Создайте модель с блоком, представляющим схему RC

Сначала создайте новую модель Simulink.

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

Использование matlabFunctionBlock для преобразования символьного результата для выхода напряжения в блок Simulink, представляющий модели объекта управления RC. matlabFunctionBlock добавляет этот новый блок в модель.

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

Добавьте больше блоков

Добавьте и расположите другие блоки с положениями и размерами относительно блока RC.

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 и v20 являются параметрами, которые нужно оценить. Во-первых, введите и инициализируйте эти параметры в рабочем пространстве MATLAB, используя начальные значения 460 μF и 1 В, соответственно. Затем создайте постоянные блоки для обоих параметров.

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);

Добавьте пандус для времени.

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

Добавьте выхода порт.

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);

Теперь блоки проводов в модели. Модель готова к 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);

Оценка параметров

Получите измеренные данные.

load sdoRCCircuit_ExperimentData

Переменные time и data загружаются в рабочую область. data - измеренное напряжение конденсатора в разы time.

Создайте sdo.Experiment объект для хранения экспериментальных данных о напряжении.

Exp = sdo.Experiment(myModel);

Создайте объект, чтобы сохранить измеренное напряжение конденсатора выхода.

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

Добавьте измеренные данные конденсатора к эксперименту в качестве ожидаемых выходных данных.

Exp.OutputData = Voltage;

Получите параметры. Установите минимальное значение для C1. Обратите внимание, что начальные предположения уже заданы.

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

Определите целевую функцию для оценки. Код для sdoRCSymbolic_Objective используемый в этом примере приведен в разделе Helper Functions в конце примера.

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

Соберите параметры модели, которые будут оценены.

v = [c1param;v20param];

Поскольку модель полностью алгебраическая, отключите предупреждающие сообщения, которые предписывают вам выбрать дискретный решатель.

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

Оцените параметры.

opt = sdo.OptimizeOptions;
opt.Method = 'lsqnonlin';
vOpt = sdo.optimize(estFcn,v,opt);
 Optimization started 25-Apr-2021 12:31:40

                                          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.

Показать предполагаемые значения.

fprintf('C1 = %e v20 = %e\n',vOpt(1).Value, vOpt(2).Value);
C1 = 2.261442e-04 v20 = 2.359446e+00

Сравнение моделируемых и экспериментальных данных

Обновите эксперимент с оцененными значениями начальных напряжений емкости и конденсатора.

Exp = setEstimatedValues(Exp,vOpt);

Симулируйте модель с оцененными значениями параметров и сравните моделируемый выход с экспериментальными данными.

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')

Figure contains an axes. The axes with title Simulated and Measured Responses contains 2 objects of type line. These objects represent Measured Voltage, Simulated Voltage.

close_system(myModel,0);

Вспомогательные функции

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