有伙伴在巴山学长交流群中询问有关如何在matlab中实现共享参数拟合的问题，感觉这个问题挺有意思的，故拿出来与大家分享。咱也根据伙伴的提问在网上进行了相关搜索，发现这个共享参数拟合的问题基本上都跟国产拟合优化神器1stopt这款软件有关。不巧咱手里正好就有1stopt的正版授权，所以就把MATLAB与1stopt两个工具的结果做一个比较。
根据1stopt帮助文件中有关共享拟合模式的概念，共享参数指的是不同的拟合公式中拥有一个或多个相同的参数。这里咱也直接从1stopt共享拟合模式的案例中copy一个简单案例，案例中自变量相同，因变量不同，公式间存在三个相同的参数，如下图所示：
上式中，m1、m2与v 为共享参数，p1、p2、p3、p4、c1、c2、c3与c4为各自独立参数。
首先给出1stopt的拟合脚本和拟合结果，脚本中SharedModel函数表明了将采用共享拟合模式来进行函数拟合：
ConstStr k1=m1*(1+1/m2), k2=m1*(1+3/m2), k3=m1*(1+10/m2), k4=m1*(1+30/m2);
SharedModel;
Variable x,y(4); //y1, y2, y3, y4;
Function y1 = v*x/(k1+x)+p1*x+c1;
y2 = v*x/(k2+x)+p2*x+c2;
y3 = v*x/(k3+x)+p3*x+c3;
y4 = v*x/(k4+x)+p4*x+c4;
Data;
//x, y1, y2, y3, y4
0. 0.0 0.0 0.0 0.0
50. 17.0261 7.227563 3.574113 0.0
100. 22.16059 13.34309 8.142564 0.1
150. 30.64281 17.76278 13.55202 4.805006
200. 33.57431 25.08648 10.82765 5.232621
400. 50.40222 38.048 21.93352 14.57796
600. 58.35754 45.00776 25.49771 13.45863
800. 62.68015 50.51803 31.82192 15.86972
1000. 63.53971 52.56842 38.12983 20.90453
10000. 81.0 74.3 70.9 61.7
1stopt拟合结果：
模型公式: y1 = v*x/(m1*(1+1/m2)+x)+p1*x+c1
y1 = 79.285487953062*x/(161.92240019*(1+1/1.69035902588378)+x)+0.000243188152287981*x+1.33059469706138
模型公式: y2 = v*x/(m1*(1+3/m2)+x)+p2*x+c2
y2 = 79.285487953062*x/(161.92240019*(1+3/1.69035902588378)+x)+(-0.000109829649274472)*x+(-0.534912686375696)
模型公式: y3 = v*x/(m1*(1+10/m2)+x)+p3*x+c3
y3 = 79.285487953062*x/(161.92240019*(1+10/1.69035902588378)+x)+(-9.86675078550804e-5)*x+0.41551196660862
模型公式: y4 = v*x/(m1*(1+30/m2)+x)+p4*x+c4
y4 = 79.285487953062*x/(161.92240019*(1+30/1.69035902588378)+x)+6.26503046551593e-5*x+0.404808726452002
迭代数: 46
计算用时(时:分:秒:微秒): 00:00:09:55
优化算法: 通用全局优化算法(UGO1)
计算结束原因: 达到收敛判断标准
均方差(RMSE): 1.55565198313939
残差平方和(SSE): 96.8021237058205
相关系数(R): 0.996911737530229
相关系数之平方(R^2): 0.99383301242554
修正R平方(Adj. R^2): 0.986133009531572
确定系数(DC): 0.993760361119802
F统计(F-Statistic): -24.9047633355909
参数 最佳估算
-------------------- -------------
v 79.285487953062
m1 161.92240019
m2 1.69035902588378
p1 0.000243188152287981
c1 1.33059469706138
p2 -0.000109829649274472
c2 -0.534912686375696
p3 -9.86675078550804E-5
c3 0.41551196660862
p4 6.26503046551593E-5
c4 0.404808726452002
不得不惊叹1stopt的通用全局优化算法的强大，在未给任何初始值的情况下仅花了9秒的时间就完成了拟合，且y1,y2,y3与y4各自拟合值与实际值的mse仅为2.2616 , 0.85414 , 2.8022 , 3.7623，如下图：
在matlab中，并没有明确共享参数拟合的这种概念，而且所有拟合均需要提供合适初始值或范围。本文决定采用直接法——非线性拟合函数lsqcurvefit与智能优化算——粒子群算法函数particleswarm来求解。
一、直接法——非线性拟合函数lsqcurvefit
为了更快达到较为理想的效果，参考1stopt的拟合结果，将所有参数的范围限定在-300到+300之间，初始值在此范围中随机生成，相应代码如下：
clear;close all;
%原始数据：x, y1, y2, y3, y4
dat = [0. 0.0 0.0 0.0 0.0
50. 17.0261 7.227563 3.574113 0.0
100. 22.16059 13.34309 8.142564 0.1
150. 30.64281 17.76278 13.55202 4.805006
200. 33.57431 25.08648 10.82765 5.232621
400. 50.40222 38.048 21.93352 14.57796
600. 58.35754 45.00776 25.49771 13.45863
800. 62.68015 50.51803 31.82192 15.86972
1000. 63.53971 52.56842 38.12983 20.90453
10000. 81.0 74.3 70.9 61.7];
t = dat(:,1); y1= dat(:,2); y2= dat(:,3); y3= dat(:,4); y4= dat(:,5);
%{
共享参数：m1、m2、v；k1~k3
独立参数：p1, p2, p3, p4, c1, c2, c3, c4; k4~k11
中间参数：k1 = m1*(1+1/m2), k2 = m1*(1+3/m2), k3 = m1*(1+10/m2), k4 = m1*(1+30/m2);
原方程：y1 = v*x/(k1+x)+p1*x+c1;
y2 = v*x/(k2+x)+p2*x+c2;
y3 = v*x/(k3+x)+p3*x+c3;
y4 = v*x/(k4+x)+p4*x+c4;
%}
fun = @(k,t) [k(3)*t./(k(1)*(1+1/k(2))+t)+k(4)*t+k(8);...
k(3)*t./(k(1)*(1+3/k(2))+t)+k(5)*t+k(9);...
k(3)*t./(k(1)*(1+10/k(2))+t)+k(6)*t+k(10);...
k(3)*t./(k(1)*(1+30/k(2))+t)+k(7)*t+k(11)];
lb = -300*ones(1,11);
ub = 300*ones(1,11);
k0 = -300+600*rand(1,11);
tic;
[kx,resnorm,residual,exitflag,output,lambda,jacobian] = lsqcurvefit(fun, k0, t,[y1;y2;y3;y4],lb,ub,[]);
disp(['lsqcurvefit用时：',num2str(toc),'秒']);
len = length(t);
ny = fun(kx,t);
ny1 = ny(1:len);
ny2 = ny(len+1:2*len);
ny3 = ny(2*len+1:3*len);
ny4 = ny(3*len+1:end);
figure('Color','w');
tx = 1:len;
plot(tx,y1,'r*',tx,ny1,'ro-','MarkerSize',7,'LineWidth',1.8);
hold on;
plot(tx,y2,'b*',tx,ny2,'bo-','MarkerSize',7,'LineWidth',1.8);
plot(tx,y3,'k*',tx,ny3,'ko-','MarkerSize',7,'LineWidth',1.8);
plot(tx,y4,'m*',tx,ny4,'mo-','MarkerSize',7,'LineWidth',1.8);
hold off;
legend('y1真实值','y1拟合值','y2真实值','y2拟合值','y3真实值','y3拟合值','y4真实值','y4拟合值','Location','best');
title(['MSE: y1: ',num2str(mse(y1,ny1)),' , y2: ',num2str(mse(y2,ny2)),' ,y3: ',num2str(mse(y3,ny3)),' ,y4: ',num2str(mse(y4,ny4))]);
非线性拟合函数lsqcurvefit的拟合结果如下图所示：
不难看出，非线性拟合函数的结果与1stopt还是存在很大差距的。需要指出的是，由于初始值是随机给定的，并不能保证每次拟合结果都完全相同，特别糟糕的情况下甚至会出现拟合失败的现象。
二、智能优化算——粒子群算法函数particleswarm
对于粒子群算法，无需给定初始，仅需给出参数的范围即可，同样采用-300到+300的区间，相关代码如下：
tic;
nvars = 11;
lb = -300*ones(1,11);
ub = 300*ones(1,11);
options = optimoptions('particleswarm','SwarmSize',100,...
'MaxIterations',200,...
'HybridFcn',@fmincon,'PlotFcn','pswplotbestf');
[kx,fval,exitflag,output,points] = particleswarm(@psofcn,nvars,lb,ub,options);
disp(['PSO用时：',num2str(toc),'秒']);
len = length(t);
ny = fun(kx,t);
ny1 = ny(1:len);
ny2 = ny(len+1:2*len);
ny3 = ny(2*len+1:3*len);
ny4 = ny(3*len+1:end);
figure('Color','w');
tx = 1:len;
plot(tx,y1,'r*',tx,ny1,'ro-','MarkerSize',7,'LineWidth',1.8);
hold on;
plot(tx,y2,'b*',tx,ny2,'bo-','MarkerSize',7,'LineWidth',1.8);
plot(tx,y3,'k*',tx,ny3,'ko-','MarkerSize',7,'LineWidth',1.8);
plot(tx,y4,'m*',tx,ny4,'mo-','MarkerSize',7,'LineWidth',1.8);
hold off;
legend('y1真实值','y1拟合值','y2真实值','y2拟合值','y3真实值','y3拟合值','y4真实值','y4拟合值','Location','best');
title(['MSE: y1: ',num2str(mse(y1,ny1)),' , y2: ',num2str(mse(y2,ny2)),' ,y3: ',num2str(mse(y3,ny3)),' ,y4: ',num2str(mse(y4,ny4))]);
function re = psofcn(x)
% 目标函数，此处选择实际值与拟合值的mse的最大值作为优化目标
dat = [0. 0.0 0.0 0.0 0.0
50. 17.0261 7.227563 3.574113 0.0
100. 22.16059 13.34309 8.142564 0.1
150. 30.64281 17.76278 13.55202 4.805006
200. 33.57431 25.08648 10.82765 5.232621
400. 50.40222 38.048 21.93352 14.57796
600. 58.35754 45.00776 25.49771 13.45863
800. 62.68015 50.51803 31.82192 15.86972
1000. 63.53971 52.56842 38.12983 20.90453
10000. 81.0 74.3 70.9 61.7];
t = dat(:,1); y1= dat(:,2); y2= dat(:,3); y3= dat(:,4); y4= dat(:,5);
m1 = x(1);
m2 = x(2);
v = x(3);
p1 = x(4);
p2 = x(5);
p3 = x(6);
p4 = x(7);
c1 = x(8);
c2 = x(9);
c3 = x(10);
c4 = x(11);
k1 = m1*(1+1/m2);
k2 = m1*(1+3/m2);
k3 = m1*(1+10/m2);
k4 = m1*(1+30/m2);
ny1 = v*t./(k1+t)+p1*t+c1;
ny2 = v*t./(k2+t)+p2*t+c2;
ny3 = v*t./(k3+t)+p3*t+c3;
ny4 = v*t./(k4+t)+p4*t+c4;
re = max([mse(y1,ny1),mse(y2,ny2),mse(y3,ny3),mse(y4,ny4)]);
end
粒子群算法函数particleswarm的拟合结果如下图所示：
相较于非线性拟合函数，粒子群优化算法具有更强的全局优化能力，结果也更接近1stopt的拟合结果。当然正如前面所讲，由于初始值都是随机给的，并不能保证每次结果都一定收敛，可能有时候效果好有时候效果差，针对这种情况，建议就是多跑几次程序。
本文案例是共享参数拟合最简单的案例，即自变量是相同的。对于自变量取值不相同、自变量长度不一、甚至是缺省等复杂情况用matlab又该如何处理呢？欢迎感兴趣的伙伴留言分享自己的见解看法。
点此阅读原文
![[工具]GitHub + PicGo 搭建免费博客图床](https://img-blog.csdnimg.cn/img_convert/8be8b771d296d9064edf3526fcce5c0c.png)
