程序请教（附EES源码）

aleisia

我要做的是通过最小二乘法，来从一组非线性方程组中确定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[i]=lookup('Lookup 1',i,'v_1')  "0.0907483"   " m^3/s"
v_2[i]=lookup('Lookup 1',i,'v_2')  "0.183561"  "m^3/s"
v_3[i]=lookup('Lookup 1',i,'v_3')  "0.715177"  "m^3/s"
v_4[i]=lookup('Lookup 1',i,'v_4')  "0.192281"  "m^3/s"
v_5[i]=lookup('Lookup 1',i,'v_5')  "0.401011" "m^3/s"

"--Wind pressure at the openings--"
p_1wa[i]=lookup('Lookup 1',i,'p_1wa')  "-47.5"   "Pa"
p_2wa[i]=lookup('Lookup 1',i,'p_2wa')   "40"
p_3nda[i]=lookup('Lookup 1',i,'p_3nda')  "32.5"
p_3nwa[i]=lookup('Lookup 1',i,'p_3nwa')  "17.5"
p_3swa[i]=lookup('Lookup 1',i,'p_3swa')   "-37.5"
p_4wa[i]=lookup('Lookup 1',i,'p_4wa')      "7.5"
p_5lwa[i]=lookup('Lookup 1',i,'p_5lwa')     "-22.5"
p_5rwa[i]=lookup('Lookup 1',i,'p_5rwa')     "-7.5"
p_6sda[i]=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[i]=lookup('Lookup 1',i,'p_1a')    "30"
p_2a[i]=lookup('Lookup 1',i,'p_2a')   "30"
p_3sa[i]=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[i]=lookup('Lookup 1',i,'p_4a')  "-22.5"
p_5a[i]=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`[i]=lookup('Lookup 1',i,'deltap_1a')    "51.80"
deltap_2a`[i]=lookup('Lookup 1',i,'deltap_2a')    "51.60"
deltap_4a`[i]=lookup('Lookup 1',i,'deltap_4a')    "53.14"
deltap_5a`[i]=lookup('Lookup 1',i,'deltap_5a')    "49.42"
deltap_6a`[i]=lookup('Lookup 1',i,'deltap_6a')    "38.12"
deltap_3ndagauge`[i]=lookup('Lookup 1',i,'deltap_3ndagauge')     "38.12"
deltap_3nwagauge`[i]=lookup('Lookup 1',i,'deltap_3nwagauge')    "38.12"
deltap_3swagauge`[i]=lookup('Lookup 1',i,'deltap_3swagauge')    "38.12"



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

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

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

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

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

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

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

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


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

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

A_3na[i]=A_6a[i]

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

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


END


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

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