6.2 有边界的线性拟合
在积累了一些进料和出料的样本数据后,我们希望计算出收率矩阵。
如果没有任何约束条件,用最小二乘法就可以得到收率W。
W=linefit(X,Y)
其中linefit(…)是最小二乘线性拟合函数。
这时的W没有任何限制,其中的元素可能大于1或小于0,这对于生产就没有意义了,所以当考虑约束1时,简单的用最小二乘法来做就行不通了。
如果只考虑第一个约束条件,0和1是wij的边界,在这个边界范围内使得误差平方和MSE最小,即:
MSE=(Y-XW)2
wij∈[0,1]
可以证明,最优解要么出现在边界上,要么出现在最小二乘法拟合的值上,也就是说wij有三种可能的取值,分别是:
wij=0
wij=1
wij=linefit(…),且满足wij∈[0,1]
当进料数有限时,穷举每种进料收率的3种可能就能拟合出最优解,m种进料就有3m种可能,找到使MSE最小的一套收率即为边界范围内的最优解。我们把这种边界约束条件下求最优解的方法称为有边界的线性拟合算法。
以第一列出料Yc1为例介绍有边界的线性拟合算法,其他进料是一样的。
1. 穷举进料收率所有可能的取值情况
每种进料的收率有3种可能取值,m种进料就相当于m重循环,所以可能的收率组合共有3m种,记为s。
s=3m
比如所有进料收率都取最小二乘拟合值时:
W(t) =linefit(X,Yc1)
其中W(t)是Yc1第t个收率矩阵,t∈[1,s]。
当所有进料收率都取下边界时:
此时的收率矩阵就是m个0组成的m*1矩阵,其他取到边界的情况是一样的,边界本身就是该进料的收率。
当只有部分进料收率取到边界时:
其他进料的收率需要最小二乘法来求,假设进料Xca,Xcb,…,Xcl的收率取到了边界,这些进料的收率就是边界,为Wra,Wrb,…,Wrl。
收率没取到边界的进料记为X’,对应的出料需要减去边界收率与进料的乘积,记为Yc1’
Yc1’= Yc1-(Xca*Wra+Xcb* Wrb+…+Xcl*Wrl)
X’的收率矩阵用最小二乘法来拟合,记为W(t)’
W(t)’=linefit(X’,Yc1’)
W(t)’与Waj,Wbj,…,Wkj共同构成第t个收率矩阵W(t)。
穷举完全部情况,共得到s个收率矩阵,从这s个收率矩阵中找出满足边界约束(即wij∈[0,1])且MSE最小的收率就是最优收率。
2. 筛选满足边界条件的收率矩阵集合[W]’
把s套收率矩阵的集合记为[W]
[W]=[W(1),W(2),…,W(s)]
[W]’=[W].select(0≤W(t)i≤1)
其中W(t)i是W(t)中第i个进料的收率。
3. 计算[W]’中每套收率与Yc1的MSE。
MSE=(Yc1-X*W(t))2
使MSE最小的W(t)即为最佳收率矩阵Wr1。
SPL例程
A |
B |
C |
D |
|
1 |
[[1,2],[1,3],[1,4]] |
/X |
||
2 |
[[0.3],[2.4],[1.6]] |
/Y |
||
3 |
[[0,1],[0,1]] |
/BD |
||
4 |
[0,1,2] |
/0:最小二乘拟合值,1:下边界,2:上边界 |
||
5 |
=(A1.~.len()-1).iterate(A4.conj((idx=#,~~.([A4(idx)].insert(0,~)))),A4) |
/3m种取值标记 |
||
6 |
=A5.(func(A12,A1,A2,A3,~)) |
/[W] |
||
7 |
=A6.select((w=~.conj(),w.id(between(~,A3(#)(1):A3(#)(2)))==[true])) |
/满足边界约束的[W]' |
||
8 |
=A7.(mse(mul(A1,~),A2)) |
/MSE |
||
9 |
=A8.pmin() |
/最小MSE的W(t)索引 |
||
10 |
=A7(A9) |
/最优W |
||
11 |
/求W(t),参数:X,Y,边界,1种取值标记 |
|||
12 |
func |
|||
13 |
=D12.align@ap([true,false],~!=0) |
|||
14 |
=B13(1) |
/第几个进料取边界 |
||
15 |
=D12(B14) |
/上界,下界,最小二乘拟合值? |
||
16 |
=B13(2) |
/第几个进料不取边界 |
||
17 |
=C12(B14) |
|||
18 |
=B17.([~(B15(#))]) |
/[Wra,Wrb,…,Wrl] |
||
19 |
if B16.len()==0 |
return B18 |
/都是边界取边界 |
|
20 |
=A12.(~(B14)) |
/取边界的[Xca,Xcb,…,Xcl] |
||
21 |
=if(B18.len()<1,B12.([null]),mul(B20,B18)) |
/Xca*Wra+Xcb* Wrb+…+Xcl*Wrl |
||
22 |
=B12.(~--B21(#)) |
/Ycj' |
||
23 |
=A12.(~(B16)) |
/不取边界的X' |
||
24 |
=w=linefit(B23,B22),if(ifnumber(w),[[w]],w) |
/W(t)' |
||
25 |
=B18|B24 |
/合并[Wra,Wrb,…,Wrl]与拟合W(t)' |
||
26 |
=B13.conj() |
/进料索引 |
||
27 |
=B26.psort() |
/进料原序 |
||
28 |
=B25(B27) |
/整理后的W(t) |
||
29 |
return B28 |
本例程介绍单出料的情况,多出料是一样的,只是多循环几次。
A5格利用迭代函数iterate(…)构造3m种可能的取值标记。
A4中的标记表示3种可能情况,0表示需要最小二乘拟合,1表示下界,2表示上界。本例中有2列进料,会有9种标记组合,如下表:
序号 |
标记组合 |
第 1 列进料 |
第 2 列进料 |
1 |
[0,0] |
最小二乘拟合 |
最小二乘拟合 |
2 |
[0,1] |
最小二乘拟合 |
下界 |
3 |
[0,2] |
最小二乘拟合 |
上界 |
4 |
[1,0] |
下界 |
最小二乘拟合 |
5 |
[1,1] |
下界 |
下界 |
6 |
[1,2] |
下界 |
上界 |
7 |
[2,0] |
上界 |
最小二乘拟合 |
8 |
[2,1] |
上界 |
下界 |
9 |
[2,2] |
上界 |
上界 |
A12格代码块计算出料Y的最佳收率矩阵W
计算结果示例:
先来看最小二乘法对上述X、Y的拟合结果:
该收率W的MSE:
MSE=0.467
再来看有边界线性拟合算法的拟合结果:
该收率W的MSE:
MSE=0.486
单纯比较MSE,最小二乘法的拟合结果更好,但因为不满足边界约束条件,不能应用在实际生产活动中,还是要选用有边界的线性拟合算法。