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

В этом примере используются:

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

В этом примере возможности символьного математического 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 = 
 $v_{1} - e^{- \frac{t}{C_{1} R_{1}}} (v_{1} - v_{20})$

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

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

v2sol = 
 $e^{- \frac{0.0001 t}{C_{1}}} (v_{20} - 5.0) + 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 $мкФ$ и 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.

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 используется в этом примере в разделе Вспомогательные функции в конце примера.

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

Документация