参数估计

计算科学 matlab 非线性规划 回归 非线性方程 最小二乘
2021-12-14 21:01:30

我必须估计一个参数(K),但我不知道我该怎么做。我认为是回归模型(最小二乘法?),但我不确定。该系统是:

dX1/dt = f(t)*X2*X1
dX2/dt = -K*f(t)*X2*X1

在哪里:

X1 和 X2 是状态变量,f(t) 是一个时变函数,K 是我要估计的参数。我知道每个采样时间 Ts=0.5 秒的 X1 和 X2 的值。我正在使用 Matlab。

4个回答

所以,有两种方法可以解决这个问题。

解决此问题的简单、不严格的方法是创建一个调用 MATLAB ODE 求解器的函数,使用ķ作为参数,并返回解决方案(X1(),X2())对应于测量数据点的所有时间。然后使用这个函数构造一个平方和误差目标函数,并使用 MATLAB 的无约束优化求解器最小化。

解决这个问题的严格方法是查看解决嵌入式动力系统优化问题的文献。据我所知,不存在将嵌入式动态系统的优化问题作为输入并返回严格全局的 MATLAB 求解器ε-最佳。此类方法通常需要非凸全局优化求解器(据我所知,MATLAB 没有)、区间算术库(MATLAB 确实存在;INTLAB 就是这样一个库)、自动微分工具(使用运算符重载的 MATLAB 也存在和面向对象的语言特性),以及大量的数学理论。

由于人们往往不想理会所有的数学理论(有很多),第一种方法往往受到从业者的青睐(这是我在本科时所学的),即使它只产生一个本地由于问题的非凸性而最优(在大多数情况下,无论如何)。有时,第二种方法对参数拟合很有用;我的两位顾问都将这种方法用于化学动力学应用,因为第一种方法是不够的。(不过,我正在努力记住论文的名称,并且时间紧迫,我稍后将不得不回到这个答案。)

卡尔曼滤波器(及其变体)可以用来解决这个问题。

你有一个所谓的灰盒系统识别问题。如果您安装了系统识别工具箱,您基本上就完成了。

http://www.mathworks.se/products/sysid/description5.html

http://www.mathworks.se/help/ident/ug/estimating-nonlinear-grey-box-models.html

如果要找到的参数是单个数字,那么您可以根据大量 K 进行蛮力预测计算(正如 Geoff 建议的那样),然后绘制计算预测和测量数据之间的差异。

这给出了改变参数将如何影响所有测量的“完整画面”,并且成本并不高,因为在现代计算机上解决数千个 ODE 只需要几分钟。然后,您可以直观地检查结果以查看要选择的 K 以及结果对 K 变化的敏感程度 - 即,您应该对结果有多确定/不确定。

此时,您可以根据需要对结果进行一些简单的最小二乘或统计分析,但您可能会通过查看结果更好地了解正在发生的事情。

这是一些进行蛮力计算和绘图的代码!我将此代码发布到公共领域。

% Some randomish input data (replace with your times tt and function f)
tt=0:.01:1;
Ktrue = 2.31; 
ff =  ifft(fft(exp(randn(1,length(tt)))).* ...
              fft(exp(-(tt-5).^2/(1.5)^2)))';
ff = ff/mean(ff);
f = @(t) interp1(tt,ff,t);

% manufacturing state data X1,X2 (replace with your X1,X2)
X1initial = 0.93; 
X2initial = 1.87;
odefuntrue = @(t,Y) [f(t)*Y(1)*Y(2), -Ktrue*f(t)*Y(1)*Y(2)]';
[T,Ytrue] = ode45(odefuntrue,tt,[X1initial,X2initial]);  
X1 = Ytrue(:,1); 
X2 = Ytrue(:,2);

% Solve the ODE for a lot of different values of K
K=0.01:0.01:5;
Ts = zeros(length(tt),length(K)); 
X1s = Ts; X2s = Ts;
for jj=1:length(K)
    odefun = @(t,Y) [f(t)*Y(1)*Y(2), -K(jj)*f(t)*Y(1)*Y(2)]';
    [T,Y] = ode45(odefun,tt,[X1(1),X2(1)]); 
    X1s(:,jj)=Y(:,1);
    X2s(:,jj)=Y(:,2);
end

% Compute the difference between the measured state, (X1,X2), and
% the predicted state for all different K's, (X1s,X2s)
errors1 = X1*ones(1,length(K)) - X1s;
errors2 = X2*ones(1,length(K)) - X2s;

% Plot the full errors at all times for different values of K
figure;
subplot(1,2,1); 
imagesc(K,tt,errors1);
title('error in X1');
xlabel('K');
ylabel('t');
colorbar;

subplot(1,2,2); 
imagesc(K,tt,errors2);
title('error in X2');
xlabel('K');
ylabel('t');
colorbar;

这是结果。您可以看到沿垂直线的误差如何为零ķ=2.31对应于真ķ, 和水平线=0,并且变为非零为ķ偏离真实值,时间增加。

在此处输入图像描述

在此处输入图像描述