6.4 质量守恒约束下的线性拟合

 

有边界的线性拟合方法保证拟合时收率被限制在边界内,误差限制法保证收率不偏离基础收率太远,现在只有约束2(所有出料对某一种进料的收率和等于1)还没有满足,本节将介绍一种线性变换的方法来满足此约束,从而求得最终的收率W

利用进料矩阵X和基础收率B求得基础出料矩阵,用YB表示。

YB=X*B

对比YYB的差距,得到两者的偏差矩阵,用E表示。

E=Y-YB

用误差限制法计算偏差收率边界BD

BD=bd(X,Y,B,rg)

其中bd(…)是误差限制法,rg是边界调节比例rg[0,1]

用有边界的线性拟合算法拟合XE,得到偏差收率,用WE表示。

WE=BDfit(X,E,BD)

其中BDfit(…)是有边界线性拟合算法。

WE不能保证W满足约束2,但可以经过线性变换,在损失准确性最小的情况下得到满足约束条件的最佳偏差收率WD

WD=L(WE)

其中L(…)是某种线性变换。

我们期望WD与基础收率B相加得到收率W,对W约束就转换为对WD的约束,如下图:

..

dijuij分别是上界和下界。

理想情况下,WE各行的总和WRL0,即WE的每一列与其他列的加和是线性关系。

利用这一线性关系可以对WE进行线性变换得到满足约束2WD

假设WE矩阵如下:

..

(1) 计算WE的各行总和WRL

..

(2) 计算WEj列以外的其他向量的总和组成的矩阵WS

..

(3) 线性变换得到WD

WEcjWScjWRL三者有如下关系:

WEcj+WScj=WRL

WEcj是有边界线性拟合算法得到的第j列最佳收率,则WScj=WRL-WEcj也是准确的,进行如下线性变换得到WD’

WD’cj=WEcj*(n-1)/n+ WScj*(-1/n)

=WEcj*(n-1)/n+(WEcj-WRL)/n

= WEcj-WRL/n

其中WD’cj是线性变换后的矩阵WD’的第j列元素。

线性变换后,相当于某个偏差收率weij减去其所在行的收率均值WRLi

..

wd’ij=weij-WRLi/n

WD’满足约束2,但仍需验证WD’能否使W满足约束1,即wij=wd’ij+bij[0,1]

如果W满足约束1,则拟合结束。

WD=WD’

如果存在wd’ij使得wij[0,1]时,需要调整变换规则。

假设WE中有z列存在不满足约束的元素。

WE=[WEc1,WEc2,…,WEc(1),…,WEc(z),…,WEcn]

其中[WEc(1),…,WEc(z)]是存在转换后仍不符合约束条件偏差收率列,即存在:

wei(p)+ bi(p)[0,1]p[1,z]

刚才的线性变换是WD’cj=WEcj*(n-1)/n+WScj*(-1/n),两者的变换系数是确定的,分别是(n-1)/n(-1/n),但因为无法满足约束,所以不能用,可以换一个思路,把WEcj的变换系数固定为1WScj的系数用解方程的方式求出来。

WScj的系数向量是向量A,方程如下:

..

其中WDcj是偏差系数矩阵WD的第j列元素。

方程①②即可求得系数向量A,还要考虑基础收率bij等于0时,对应的收率weij也要等于0

解方程得到向量A

..

有了向量A后,带入方程①中,求得WDcj

此时还要再验证一遍WD是否能使W满足约束1,如果满足,则计算结束,求得WD。否则还要把不满足约束条件的那一列加入到[WEc(1),…,WEc(z)],重复上述方法,直到所有wij都满足约束1wij=(wdij+bij)[0,1]),得到WD

最后用WD加基础收率矩阵B得到W

W=B+WD

SPL例程


A

B

C

D

1

[[30,8],[31,7],[38,10]]

/X


2

[[2,13,23],[3,15,20],[11,13,24]]

/Y


3

[[0,0.5,0.5],[0.55,0.05,0.4]]

/B


4

0.1

/rg


5

=mul(A1,A3)

/X*B


6

=A2.(~--A5(#))

/E


7

=bd(A1,A2,A3,A4).(~.(~.(round(~,3))))

/BD


8

=A6.~.((idx=#,bdfit(A1,A6.([~(idx)]),A7.(~(idx))).conj()))

/WET

9

=transpose(A8)

/WE


10

=func(A12,A9,A3)

/WD


11

/线性变换,参数:WEB



12

func




13


=A12.(~.count(~==0))


/WE每行0的个数

14


=A12.~.len()


/出料数

15


=A12.((s=~.sum(),ind=#,~.(if(~==0,0,s/(B14-B13(ind))))))

/各行求均值

16


=A12.(~--B15(#))

/weij-WRLi/n

17


=B12.(~++B16(#))

/W=WD'+B

18


=transpose(B17)


19


=B18.pselect@a(~.pselect(~<0||~>1)>0)

/记录WD‘不满足约束的列索引

20


if B19.len()==0

return B17

/没有超限返回W

21


else

=func(A24,A12,B12,B19)

/超限时变换规则

22



return C21

23

/递归消除不符合约束2的元素

24

func




25


=B24.~.len()


26


=A24.(~.sum())


/WRL

27


=A24.((w=~,w.((ind=#,if(~==0,0,w(to(B25).delete(ind)).sum())))))

/WS

28


=A24.(~.pselect@a(~==0))

/WEri中收率为0的索引

29


=B28.(~&C24)

/不满足约束及等于0的索引

30


=B26.(~/((B25-B29(#).len()-1)*~+A24(#)(B29(#)).sum()))

/ai

31


=A24.((ind=#,~.(if(B29(ind).pos(#)>0,0,B30(ind)))))

/A

32


=B27.(~**B31(#))


/A**WScj

33


=A24.((idx=#,if(B31(#).id()==[0],B32(#),~--B32(#))))

/WD

34


=B24.(~++B33(#))


/B+WD

35


=transpose(B34)


/WT

36


=B35.pselect@a(~.pselect(~<0||~>1)>0)

/再次确认是否有超限收率

37


if B36.len()==0

return B34

/没有超限返回W

38


else

=C24=C24&B36

/有超限记录下来

39


=func(A24,A24,B12,C24)


/WD

A7格中的bd(…)是误差限制法,用来计算偏差收率边界BD

A8格中的bdfit(…)是有边界的线性拟合算法。

A12格代码块对WE做线性变换;

A24格代码块对WD’做变换;

计算结果示例:

进料X

..

出料Y

..

基础收率B

..

边界调节比例rg

rg=0.1

算出的边界BD

..

最终的收率W

..

先来看使用基础收率B,各出料的MSE

MSE1=12.24

MSE2=16.24

MSE3=8.98

其中MSEj是第j个出料的MSE

再看使用W,各出料的MSE

MSE1=10.97

MSE2=5.13

MSE3=3.86

很明显是拟合后的W效果更好。