非线性方程组约束下的最小二乘法拟合问题

aleisia

大家好！
我刚开始学习1stopt，主要参考的就是官方的教程，可是教程真的很简单，有些地方不是很确定的感觉。

我要做的是通过最小二乘法，来从一组非线性方程组中确定5个常数。
具体的目标函数就是下面程序里的SSRE_total,它的含义就是实验值和理论值之间的平方差，我有98组数据，所以要从一个数据表中读入一些实验值。
有的压力是实验测出来的，再通过下面的程序算出来这些压力（计算值），具体计算就是基于连续性方程和假定的一个压力与流量，面积的关系式v=c*deltaP^n*A.

不知道高手能不能帮忙写成1stopt的程序？

我用Engineering Equation Solver写的程序如下：



c_2=0.55

c_out=180.99    "mg/m^3"
eta=0.999975
c_f=(1-eta)*c_out

{c=0.45
n=0.7
d_0=0.02653
d_1=0.01909
d_2=0.001306}

M=98
DUPLICATE i=1,M                        

v_1=lookup('Lookup 1',i,'v_1')  "0.0907483"   " m^3/s"
v_2=lookup('Lookup 1',i,'v_2')  "0.183561"  "m^3/s"
v_3=lookup('Lookup 1',i,'v_3')  "0.715177"  "m^3/s"
v_4=lookup('Lookup 1',i,'v_4')  "0.192281"  "m^3/s"
v_5=lookup('Lookup 1',i,'v_5')  "0.401011" "m^3/s"

"--Wind pressure at the openings--"
p_1wa=lookup('Lookup 1',i,'p_1wa')  "-47.5"   "Pa"
p_2wa=lookup('Lookup 1',i,'p_2wa')   "40"
p_3nda=lookup('Lookup 1',i,'p_3nda')  "32.5"
p_3nwa=lookup('Lookup 1',i,'p_3nwa')  "17.5"
p_3swa=lookup('Lookup 1',i,'p_3swa')   "-37.5"
p_4wa=lookup('Lookup 1',i,'p_4wa')      "7.5"
p_5lwa=lookup('Lookup 1',i,'p_5lwa')     "-22.5"
p_5rwa=lookup('Lookup 1',i,'p_5rwa')     "-7.5"
p_6sda=lookup('Lookup 1',i,'p_6sda')     "-35"

"--Wind pressure at the position of pressure transducer--"
"For zone 3, the two windward measurement positions are the same as the above openings, so the wind pressure values are the same"
p_1a=lookup('Lookup 1',i,'p_1a')    "30"
p_2a=lookup('Lookup 1',i,'p_2a')   "30"
p_3sa=lookup('Lookup 1',i,'p_3sa')    "-35"    "p_3sa is not the same as p_3swa, p_3sa is approximated the same as p_6sda, because the two transducers are close to each other, but a little far from south window in zone 3"
p_4a=lookup('Lookup 1',i,'p_4a')  "-22.5"
p_5a=lookup('Lookup 1',i,'p_5a')  "-22.5"

"--Pressure difference measured at the position of the pressure transducer, not necessarily at the openings--"
deltap_1a`=lookup('Lookup 1',i,'deltap_1a')    "51.80"
deltap_2a`=lookup('Lookup 1',i,'deltap_2a')    "51.60"
deltap_4a`=lookup('Lookup 1',i,'deltap_4a')    "53.14"
deltap_5a`=lookup('Lookup 1',i,'deltap_5a')    "49.42"
deltap_6a`=lookup('Lookup 1',i,'deltap_6a')    "38.12"
deltap_3ndagauge`=lookup('Lookup 1',i,'deltap_3ndagauge')     "38.12"
deltap_3nwagauge`=lookup('Lookup 1',i,'deltap_3nwagauge')    "38.12"
deltap_3swagauge`=lookup('Lookup 1',i,'deltap_3swagauge')    "38.12"



"--Absolute pressure in each zone--"
p_1=deltap_1a+p_1a
p_2=deltap_2a+p_2a

{p_3_ave=(p_3nda+deltap_3nda+p_3nwa+deltap_3nwa+p_3sa+deltap_3swa)/3         "absolute average pressure in zone 3"}
p_4=deltap_4a+p_4a
p_5=deltap_5a+p_5a
p_6=deltap_6a+p_6sda

"--Pressure difference averaged for the 3 openings in zone 3--"
deltap_3nda=IF(p_3_ave, p_3nda, p_3nda-p_3_ave, 0, p_3_ave-p_3nda)    "105.0"
deltap_3nwa=IF(p_3_ave, p_3nwa, p_3nwa-p_3_ave, 0, p_3_ave-p_3nwa)   "95.19"         "90.59=(105+95.19+71.59)/3?"
deltap_3swa=IF(p_3_ave, p_3swa, p_3swa-p_3_ave, 0, p_3_ave-p_3swa)   "71.59"

"--Pressure difference at the opening--"
deltap_1wa=IF(p_1, p_1wa, p_1wa-p_1, 0, p_1-p_1wa)
deltap_2wa=IF(p_2, p_2wa, p_2wa-p_2, 0, p_2-p_2wa)
deltap_4wa=IF(p_4, p_4wa, p_4wa-p_4, 0, p_4-p_4wa)
deltap_5lwa=IF(p_5, p_5lwa, p_5lwa-p_5, 0, p_5-p_5lwa)
deltap_5rwa=IF(p_5, p_5rwa, p_5rwa-p_5, 0, p_5-p_5rwa)

deltap_31=IF(p_3_ave, p_1, p_1-p_3_ave, 0, p_3_ave-p_1)  
deltap_32=IF(p_3_ave, p_2, p_2-p_3_ave, 0, p_3_ave-p_2)
deltap_34=IF(p_3_ave, p_4, p_4-p_3_ave, 0, p_3_ave-p_4)  
deltap_35=IF(p_3_ave, p_5, p_5-p_3_ave, 0, p_3_ave-p_5)  
deltap_36=IF(p_3_ave, p_6, p_6-p_3_ave, 0, p_3_ave-p_6)  

{(v_1a/(c*A_1a))^(1/n)=deltap_1wa
(v_2a/(c*A_2a))^(1/n)=deltap_2wa
(v_3swa/(c*A_3swa))^(1/n)=deltap_3swa
(v_3nda/(c*A_3na))^(1/n)=deltap_3nda
(v_3nwa/(c*A_3nwa))^(1/n)=deltap_3nwa
(v_4a/(c*A_4a))^(1/n)=deltap_4wa
(v_5la/(c*A_5la))^(1/n)=deltap_5lwa
(v_5ra/(c*A_5ra))^(1/n)=deltap_5rwa
(v_6a/(c*A_6a))^(1/n)=deltap_6sda}

v_1a=c*abs(deltap_1wa)^n*A_1a
v_2a=c*abs(deltap_2wa)^n*A_2a
v_3swa=c*abs(deltap_3swa)^n*A_3swa
v_3nda=c*abs(deltap_3nda)^n*A_3na
v_3nwa=c*abs(deltap_3nwa)^n*A_3nwa
v_4a=c*abs(deltap_4wa)^n*A_4a
v_5la=c*abs(deltap_5lwa)^n*A_5la
v_5ra=c*abs(deltap_5rwa)^n*A_5ra
v_6a=c*abs(deltap_6a)^n*A_6a

v_31=c*abs(deltap_31)^n*A_31
v_32=c*abs(deltap_32)^n*A_32
v_34=c*abs(deltap_34)^n*A_34
v_35=c*abs(deltap_35)^n*A_35
v_36=c*abs(deltap_36)^n*A_36


v_1a_star=IF(p_1, p_1wa, v_1a, 0, -v_1a)
v_2a_star=IF(p_2, p_2wa, v_2a, 0, -v_2a)
v_3nda_star=IF(p_3_ave, p_3nda, v_3nda, 0, -v_3nda)
v_3nwa_star=IF(p_3_ave, p_3nwa, v_3nwa, 0, -v_3nwa)
v_3swa_star=IF(p_3_ave, p_3swa, v_3swa, 0, -v_3swa)
v_4a_star=IF(p_4, p_4wa, v_4a, 0, -v_4a)
v_5la_star=IF(p_5, p_5lwa, v_5la, 0, -v_5la)
v_5ra_star=IF(p_5, p_5rwa, v_5ra, 0, -v_5ra)
v_6a_star=IF(p_6, p_6sda, v_6a, 0, -v_6a)
v_31_star=IF(p_1, p_3_ave, v_31, 0, -v_31)
v_32_star=IF(p_2, p_3_ave, v_32, 0, -v_32)
v_34_star=IF(p_4, p_3_ave, v_34, 0, -v_34)
v_35_star=IF(p_5, p_3_ave, v_35, 0, -v_35)
v_36_star=IF(p_6, p_3_ave, v_36, 0, -v_36)

{These 6 equations are airflow mass conservation for unbalanced dilution ventilation system only}
v_1a_star+v_31_star+v_1=0
v_2a_star+v_32_star+v_2=0
v_3swa_star+v_3nda_star+v_3nwa_star+v_3=v_31_star+v_32_star+v_34_star+v_35_star+v_36_star
v_4a_star+v_34_star+v_4=0
v_5la_star+v_5ra_star+v_35_star+v_5=0
v_36_star+v_6a_star=0

A_3na=A_6a

A_1a=A_2a
A_1a=A_3nwa
A_1a=A_3swa
A_1a=A_4a
A_1a=A_5la
A_1a=A_5ra
A_31=A_32
A_31=A_34
A_31=A_35
A_31=A_36

A_1a=d_0
A_3na=d_1                                    "assuming constant crack area for the building envelop openings"
A_31=d_2


END


SSRE_1a=sum((deltap_1a`-deltap_1a)^2/deltap_1a`^2,i=1,M)/M        {sum of square errors}
SSRE_2a=sum((deltap_2a`-deltap_2a)^2/deltap_2a`^2,i=1,M)/M
SSRE_4a=sum((deltap_4a`-deltap_4a)^2/deltap_4a`^2,i=1,M)/M
SSRE_5a=sum((deltap_5a`-deltap_5a)^2/deltap_5a`^2,i=1,M)/M
SSRE_6a=sum((deltap_6a`-deltap_6a)^2/deltap_6a`^2,i=1,M)/M
SSRE_3nda=sum((deltap_3ndagauge`-deltap_3nda)^2/deltap_3ndagauge`^2,i=1,M)/M
SSRE_3nwa=sum((deltap_3nwagauge`-deltap_3nwa)^2/deltap_3nwagauge`^2,i=1,M)/M
SSRE_3swa=sum((deltap_3swagauge`-deltap_3swa)^2/deltap_3swagauge`^2,i=1,M)/M

SSRE_total=SSRE_1a+SSRE_2a+SSRE_4a+SSRE_5a+SSRE_6a+SSRE_3nda+SSRE_3nwa+SSRE_3swa