参数估计时将模型设置为稳态(代码)

这个例子使用了:

打开脚本

本例展示了在参数估计过程中如何将模型设置为稳态。在电力系统和飞机动力学等许多应用中，将模型设置为稳态非常重要。本例使用种群动态模型。

本例需要Simulink®Control D金宝appesign™软件。

模型描述

Simu金宝applink模型sdoPopulationInflux建立一个简单的生态模型，其中生物种群的增长受环境承载能力的限制。

$\frac{y}{dt} = R (1 - \frac{y}{K}) (y + 10)$

是生物种群的固有生长速率。

是环境的承载能力。

也有来自邻近环境的其他生物成员的涌入。该模型使用标准化单位。

打开模型。

modelName =“sdoPopulationInflux”;open_system (modelName)

该文件sdoPopulationInflux_Data.mat包含环境中的种群数据以及来自邻近环境的其他生物的流入数据。

负载sdoPopulationInflux_Data.mat;时间序列:Population_ts Inflow_tshFig =图;次要情节(2,1,1);情节(Population_ts)次要情节(2,1,2);情节(Inflow_ts)

种群从稳定的状态开始。一段时间后，周围环境的生物大量涌入。根据测量数据，您可以估计模型参数的值。

的参数R代表有机体固有的生长速度。使用1作为该参数的初始猜测。它是非负的。

R = sdo.getParameterFromModel(modelName，“R”）;R.Value = 1;r .最小值= 0;

的参数K表示环境的承载能力。使用2作为该参数的初始猜测。据了解，它至少是0.1。

K = sdo.getParameterFromModel(modelName，“K”）;K.Value = 2;k .最小值= 0.1;

将这些参数收集到一个向量中。

v = [R K];

比较实测数据与初始模拟输出

创建一个实验对象。

Exp = sdo.Experiment(modelName);

联系Population_ts用模型输出。

Population = 金宝appSimulink.SimulationData.Signal;人口。Name =“人口”;人口。BlockPath = [modelName .' /集成商'];人口。PortType =“输出港”;人口。PortIndex = 1;人口。Values = Population_ts;

添加人口为了实验。

Exp.OutputData =人口;

联系Inflow_ts模型输入。

流入= Simul金宝appink.SimulationData.Signal;流入。Name =“人口”;流入。BlockPath = [modelName .' / in '];流入。PortType =“输出港”;流入。PortIndex = 1;流入。值= Inflow_ts;

添加流入为了实验。

Exp.InputData =流入;

利用实验创建仿真场景，并获得仿真输出。

Exp = setEstimatedValues(Exp, v);%使用参数/状态向量模拟器= createSimulator(Exp);模拟器= sim(模拟器);

在记录的仿真数据中搜索模型输出信号。

SimLog = find(模拟器。LoggedData,…get_param (modelName“SignalLoggingName”));PopulationSim = find(SimLog，“人口”）;

模型输出与数据匹配得不是很好，这表明需要对模型参数进行更好的估计。

clf;情节(PopulationSim。值,“r”）;持有在;情节(Population_ts“b”）;传奇(“仿真模型”，的测量数据，“位置”，“最佳”）;

估计参数

为了估计参数，定义一个目标函数来计算模型模拟与实测数据之间的平方和误差。

estFcn = @(v) sdoPopulationInflux_Objective(v，模拟器，Exp);类型sdoPopulationInflux_Objective.m

function vals = sdoPopulationInflux_Objective(v, Simulator, Exp, OpPointSetup) %比较模型输出与数据% %输入:% v -参数和/或状态向量%模拟器-用于模拟模型% Exp -实验对象% OpPointSetup -对象用于设置稳态计算%工作点%版权2018 the MathWorks, Inc. %如果nargin < 4 OpPointSetup = [];需求设置req = sdo.requirements.SignalTracking;要求的事情。Type =' ==';要求的事情。方法= '残差';%模拟模型Exp = setEstimatedValues(Exp, v);%使用参数/状态的矢量模拟器=创建模拟器(Exp，模拟器);strOT = mat2str(Exp.OutputData(1).Values.Time);if isempty(OpPointSetup)模拟器= sim(模拟器，'OutputOption'， 'AdditionalOutputTimes'， 'OutputTimes'， strOT); else Simulator = sim(Simulator, 'OutputOption', 'AdditionalOutputTimes', 'OutputTimes', strOT, ... 'OperatingPointSetup', OpPointSetup); end % Compare model output with data SimLog = find(Simulator.LoggedData, ... get_param(Exp.ModelName, 'SignalLoggingName') ); OutputModel = find(SimLog, 'Population'); Model_Error = evalRequirement(req, OutputModel.Values, Exp.OutputData.Values); vals.F = Model_Error;

定义优化选项。％opts = sdo.OptimizeOptions;选择。方法=“lsqnonlin”;

估计参数。

vOpt = sdo。优化(estFcn, v, opts);disp (vOpt)

一阶Iter f -count f(x)步长最优性05 12.485 1 1 10 1.09824 1.184 0.244 2 15 0.9873 1.088 0.0259 3 20 0.952948 1.217 0.00624 4 25 0.946892 0.9151 0.00197 5 30 0.946484 0.3541 0.00153局部最小值。Lsqnonlin停止了，因为最终平方和相对于其初始值的变化小于函数公差的值。(1,1) =名称:'R'值:5.5942最小值:0最大值:无穷大Free: 1刻度:1信息:[1x1 struct](1,2) =名称:'K'值:3.2061最小值:0.1000最大值:无穷大Free: 1刻度:2信息:[1x1 struct] 1x2参数。连续

利用模型中估计的参数值，得到模型响应。在记录的仿真数据中搜索模型输出信号。

Exp = setEstimatedValues(Exp, vOpt);模拟器= createsemulator (Exp，模拟器);模拟器= sim(模拟器);SimLog = find(模拟器。LoggedData,…get_param (modelName“SignalLoggingName”));PopulationSim = find(SimLog，“人口”）;

将实测的总体数据与优化后的模型响应进行比较，仍然不能很好地匹配。在模型响应的开始有一个瞬态，在那里它与实测数据有显著的不同。

clf;情节(PopulationSim。值,“r”）;持有在;情节(Population_ts“b”）;传奇(“仿真模型”，的测量数据，“位置”，“最佳”）;

在估计过程中将模型置于稳态

为了提高模型与实测数据的拟合度，在进行参数估计时，需要将模型设置为稳态。定义一个工作点规范。输入由实验数据可知。因此，(1)计算稳态工作点时不应将其作为自由变量;(2)找到工作点后，在模拟模型时不应使用其输入。另一方面，在模拟模型时应使用计算工作点时得到的所有状态。创建一个sdo。OperatingPointSetup对象来收集工作点规范、要使用的输入和要使用的状态，这样这些信息就可以传递给目标函数，并在模拟模型时使用。您还可以提供第四个参数sdo。OperatingPointSetup指定计算工作点的选项。例如，选项“graddescent-proj”通常用于为使用物理建模的系统查找操作点。

opSpec = operspec(modelName);opSpec.Inputs(1)。已知=真实;inputsToUse = [];statesToUse = 1: number (opSpec.States);OpPointSetup = sdo。operationingpointsetup (opSpec, inputsToUse, statesToUse);

对参数进行估计，使模型在过程中处于稳态。

estFcn = @(v) sdoPopulationInflux_Objective(v，模拟器，Exp, OpPointSetup);vOpt = sdo。优化(estFcn, v, opts);disp (vOpt)

优化开始26- 2月2022 18:06:04一阶Iter f -count f(x)步长优化05 11.1517 11 10 0.025674 0.5732 0.045 2 15 0.00239293 0.3451 0.357 3 20 0.000692478 0.0148 0.00301 4 25 0.00069236 6.539e-05 1.16e-07发现局部最小值。优化完成，因为梯度的大小小于最优性公差的值。(1,1) =名称:'R'值:0.5953最小值:0最大值:无穷大Free: 1刻度:1信息:[1x1 struct](1,2) =名称:'K'值:3.0988最小值:0.1000最大值:无穷大Free: 1刻度:2信息:[1x1 struct] 1x2参数。连续