带有非线性方程组约束的最小二乘法问题
本帖最后由 denny2651 于 2012-4-19 15:59 编辑大家好！
我要做的是通过最小二乘法，来从一组非线性方程组中确定5个常数(c,n,d_0,d_1,d_2)。
具体的目标函数就是下面程序里的SSRE_total,它的含义就是实验值和理论值之间的平方差，实验数据已经通过constant定义成矩阵。
有的压力是实验测出来的，再通过下面的程序算出来这些压力（计算值），具体计算就是基于连续性方程和假定的一个压力与流量，面积的关系式v=c*deltaP^n*A. A有三种值，分别就是d_0,d_1,d_2.
下面是我稍微改动以后的1stopt程序，constant定义的就是那些原始数据，每个数组有98个原始数据。
我要处理的问题很简单，一个有6个房间的一层建筑，已知其中5个房间的供气流量和与周围环境的压差读数，需要通过假设一个leakage function，以及利用质量守恒，来估算每个房间每个裂缝的等效面积。
公式呢，最关键是质量守恒：
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;
（v_1，v_2，v_3，v_4，v_5已知）
以及流过每个房间缝隙的流量与压差，面积的一个关系，就是leakage function：比如：
v_36=c*deltap_36^n*A_36，其他的类同。
这里我假定A_36=d_2
其他的还有另两种面积，分别为d_0和d_1，加上c和n，总共是5个，需要通过最小化目标函数SSRE_total得到这5个参数。
程序的其他部分，比如那些if， else if之类的，都只是为了判断一个房间两侧压力不同时，如何把压差给计算成正值而不是负值。还有那些各种压力之间的转换，就是考虑到压力测量点和裂缝开口处的外界风压不同，需要转换转换“基准值”而已，最终基准都转换到了以大气压为准。
v_36_star和v_36的区别就是压差不同决定了在上述的质量守恒里，是不是加个"-"号，这样就使得不管两侧压力大小如何，上述的质量守恒方程形式，符号，总是对的。
整个求解的思路就是如此，这里就是通过求解质量守恒，求解出一组“计算压差值”，再与“实验压差值”比较，就是那个目标函数SSRE_total，最小化以后，希望能确定c,n,d_1,d_2,d_3这5个参数。
我的想法是：把所有的那些连续性方程（质量守恒）作为conststr，要最优化的函数就是实验值和计算值的加总平方差。可是程序提示我说：d_0,d_1,d-2没出现在minfunction里，所以无法继续。
不知道大家有什么好的方法可以解决这个问题的么？谢谢！
我写的1stopt程序如下：
Title "98GA";
Parameter d_0[-1,1],d_1[-1,1],d_2[-1,1],c,n;
Var i: integer;
constant v_1(1:98)=;
constant v_2(1:98)=;
constant v_3(1:98)=;
constant v_4(1:98)=;
constant v_5(1:98)=;
constant p_1wa(1:98)=[-0.3,-0.3,-0.3,-0.3,-0.3,-0.3,-0.3,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,-0.3,-0.3,-0.3,-0.3,-0.3,-0.3,-0.3,0,0,0,0,0,0,0,0,0,0,0,0,0,0,-0.3,-0.3,-0.3,-0.3,-0.3,-0.3,-0.3,0,0,0,0,0,0,0,0,0,0,0,0,0,0,-0.3,-2.3,-2.3,-2.3,-2.3,-2.3,-2.3,-2.3,-2.3,-2.3,-2.3,-2.3,-2.3,-2.3,-2.3,-2.3,-2.3,-2.3,-2.3,-2.3,-2.3,-2.3,-0.3,-0.3,-0.3,-0.3,-0.3,-0.3];
constant p_2wa(1:98)=;
constant p_3nda(1:98)=[-0.1,-0.1,-0.1,-0.1,-0.1,-0.1,-0.1,0,0,0,0,0,0,0,-0.1,-0.1,-0.1,-0.1,-0.1,-0.1,-0.1,0,0,0,0,0,0,0,-0.1,-0.1,-0.1,-0.1,-0.1,-0.1,-0.1,-0.1,-0.1,-0.1,-0.1,-0.1,-0.1,-0.1,0,0,0,0,0,0,0,-0.1,-0.1,-0.1,-0.1,-0.1,-0.1,-0.1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,-0.1,-2.9,-2.9,-2.9,-2.9,-2.9,-2.9,-2.9,-2.9,-2.9,-2.9,-2.9,-2.9,-2.9,-2.9,-2.9,-2.9,-2.9,-2.9,-2.9,-2.9,-2.9,-0.1,-0.1,-0.1,-0.1,-0.1,-0.1];
constant p_3nwa(1:98)=[-0.8,-0.8,-0.8,-0.8,-0.8,-0.8,-0.8,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,-0.8,-0.8,-0.8,-0.8,-0.8,-0.8,-0.8,0,0,0,0,0,0,0,0,0,0,0,0,0,0,-0.8,-0.8,-0.8,-0.8,-0.8,-0.8,-0.8,0,0,0,0,0,0,0,0,0,0,0,0,0,0,-0.8,-8.4,-8.4,-8.4,-8.4,-8.4,-8.4,-8.4,-8.4,-8.4,-8.4,-8.4,-8.4,-8.4,-8.4,-8.4,-8.4,-8.4,-8.4,-8.4,-8.4,-8.4,-0.8,-0.8,-0.8,-0.8,-0.8,-0.8];
constant p_3swa(1:98)=[-0.2,-0.2,-0.2,-0.2,-0.2,-0.2,-0.2,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,-0.2,-0.2,-0.2,-0.2,-0.2,-0.2,-0.2,0,0,0,0,0,0,0,0,0,0,0,0,0,0,-0.2,-0.2,-0.2,-0.2,-0.2,-0.2,-0.2,0,0,0,0,0,0,0,0,0,0,0,0,0,0,-0.2,-2.2,-2.2,-2.2,-2.2,-2.2,-2.2,-2.2,-2.2,-2.2,-2.2,-2.2,-2.2,-2.2,-2.2,-2.2,-2.2,-2.2,-2.2,-2.2,-2.2,-2.2,-0.2,-0.2,-0.2,-0.2,-0.2,-0.2];
constant p_4wa(1:98)=[-0.2,-0.2,-0.2,-0.2,-0.2,-0.2,-0.2,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,-0.2,-0.2,-0.2,-0.2,-0.2,-0.2,-0.2,0,0,0,0,0,0,0,0,0,0,0,0,0,0,-0.2,-0.2,-0.2,-0.2,-0.2,-0.2,-0.2,0,0,0,0,0,0,0,0,0,0,0,0,0,0,-0.2,-1.1,-1.1,-1.1,-1.1,-1.1,-1.1,-1.1,-1.1,-1.1,-1.1,-1.1,-1.1,-1.1,-1.1,-1.1,-1.1,-1.1,-1.1,-1.1,-1.1,-1.1,-0.2,-0.2,-0.2,-0.2,-0.2,-0.2];
constant p_5lwa(1:98)=[-0.6,-0.6,-0.6,-0.6,-0.6,-0.6,-0.6,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,-0.6,-0.6,-0.6,-0.6,-0.6,-0.6,-0.6,0,0,0,0,0,0,0,0,0,0,0,0,0,0,-0.6,-0.6,-0.6,-0.6,-0.6,-0.6,-0.6,0,0,0,0,0,0,0,0,0,0,0,0,0,0,-0.6,-2.2,-2.2,-2.2,-2.2,-2.2,-2.2,-2.2,-2.2,-2.2,-2.2,-2.2,-2.2,-2.2,-2.2,-2.2,-2.2,-2.2,-2.2,-2.2,-2.2,-2.2,-0.6,-0.6,-0.6,-0.6,-0.6,-0.6];
constant p_5rwa(1:98)=[-0.8,-0.8,-0.8,-0.8,-0.8,-0.8,-0.8,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,-0.8,-0.8,-0.8,-0.8,-0.8,-0.8,-0.8,0,0,0,0,0,0,0,0,0,0,0,0,0,0,-0.8,-0.8,-0.8,-0.8,-0.8,-0.8,-0.8,0,0,0,0,0,0,0,0,0,0,0,0,0,0,-0.8,-1.8,-1.8,-1.8,-1.8,-1.8,-1.8,-1.8,-1.8,-1.8,-1.8,-1.8,-1.8,-1.8,-1.8,-1.8,-1.8,-1.8,-1.8,-1.8,-1.8,-1.8,-0.8,-0.8,-0.8,-0.8,-0.8,-0.8];
constant p_6sda(1:98)=[-0.3,-0.3,-0.3,-0.3,-0.3,-0.3,-0.3,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,-0.3,-0.3,-0.3,-0.3,-0.3,-0.3,-0.3,0,0,0,0,0,0,0,0,0,0,0,0,0,0,-0.3,-0.3,-0.3,-0.3,-0.3,-0.3,-0.3,0,0,0,0,0,0,0,0,0,0,0,0,0,0,-0.3,-1.5,-1.5,-1.5,-1.5,-1.5,-1.5,-1.5,-1.5,-1.5,-1.5,-1.5,-1.5,-1.5,-1.5,-1.5,-1.5,-1.5,-1.5,-1.5,-1.5,-1.5,-0.3,-0.3,-0.3,-0.3,-0.3,-0.3];
constant p_1a(1:98)=[-0.5,-0.5,-0.5,-0.5,-0.5,-0.5,-0.5,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,-0.5,-0.5,-0.5,-0.5,-0.5,-0.5,-0.5,0,0,0,0,0,0,0,0,0,0,0,0,0,0,-0.5,-0.5,-0.5,-0.5,-0.5,-0.5,-0.5,0,0,0,0,0,0,0,0,0,0,0,0,0,0,-0.5,-2,-2,-2,-2,-2,-2,-2,-2,-2,-2,-2,-2,-2,-2,-2,-2,-2,-2,-2,-2,-2,-0.5,-0.5,-0.5,-0.5,-0.5,-0.5];
constant p_2a(1:98)=[-0.9,-0.9,-0.9,-0.9,-0.9,-0.9,-0.9,0,0,0,0,0,0,0,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0,0,0,0,0,0,0,-0.9,-0.9,-0.9,-0.9,-0.9,-0.9,-0.9,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0,0,0,0,0,0,0,-0.9,-0.9,-0.9,-0.9,-0.9,-0.9,-0.9,0,0,0,0,0,0,0,0,0,0,0,0,0,0,-0.9,-1.7,-1.7,-1.7,-1.7,-1.7,-1.7,-1.7,-1.7,-1.7,-1.7,-1.7,-1.7,-1.7,-1.7,-1.7,-1.7,-1.7,-1.7,-1.7,-1.7,-1.7,-0.9,-0.9,-0.9,-0.9,-0.9,-0.9];
constant p_3sa(1:98)=[-0.2,-0.2,-0.2,-0.2,-0.2,-0.2,-0.2,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,-0.2,-0.2,-0.2,-0.2,-0.2,-0.2,-0.2,0,0,0,0,0,0,0,0,0,0,0,0,0,0,-0.2,-0.2,-0.2,-0.2,-0.2,-0.2,-0.2,0,0,0,0,0,0,0,0,0,0,0,0,0,0,-0.2,-1.9,-1.9,-1.9,-1.9,-1.9,-1.9,-1.9,-1.9,-1.9,-1.9,-1.9,-1.9,-1.9,-1.9,-1.9,-1.9,-1.9,-1.9,-1.9,-1.9,-1.9,-0.2,-0.2,-0.2,-0.2,-0.2,-0.2];
constant p_4a(1:98)=[-0.2,-0.2,-0.2,-0.2,-0.2,-0.2,-0.2,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,-0.2,-0.2,-0.2,-0.2,-0.2,-0.2,-0.2,0,0,0,0,0,0,0,0,0,0,0,0,0,0,-0.2,-0.2,-0.2,-0.2,-0.2,-0.2,-0.2,0,0,0,0,0,0,0,0,0,0,0,0,0,0,-0.2,-1.9,-1.9,-1.9,-1.9,-1.9,-1.9,-1.9,-1.9,-1.9,-1.9,-1.9,-1.9,-1.9,-1.9,-1.9,-1.9,-1.9,-1.9,-1.9,-1.9,-1.9,-0.2,-0.2,-0.2,-0.2,-0.2,-0.2];
constant p_5a(1:98)=;
constant deltap_1a_exp(1:98)=;
constant deltap_2a_exp(1:98)=;
constant deltap_4a_exp(1:98)=;
constant deltap_5a_exp(1:98)=;
constant deltap_6a_exp(1:98)=;
constant deltap_3ndagauge_exp(1:98)=;
constant deltap_3nwagauge_exp(1:98)=;
constant deltap_3swagauge_exp(1:98)=;
conststr
for i:=1 to 98 do
p_1=deltap_1a+p_1a;
p_2=deltap_2a+p_2a;
p_4=deltap_4a+p_4a;
p_5=deltap_5a+p_5a;
p_6=deltap_6a+p_6sda;
if p_3_ave<p_3nda then deltap_3nda:=p_3nda-p_3_ave;
else if p_3_ave=p_3nda then deltap_3nda:=0;
else if p_3_ave>p_3nda then deltap_3nda:=p_3_ave-p_3nda;
if p_3_ave<p_3nwa then deltap_3nwa:=p_3nwa-p_3_ave;
else if p_3_ave=p_3nwa then deltap_3nwa:=0;
else if p_3_ave>p_3nwa then deltap_3nwa:=p_3_ave-p_3nwa;
if p_3_ave<p_3swa then deltap_3swa:=p_3swa-p_3_ave;
else if p_3_ave=p_3swa then deltap_3swa:=0;
else if p_3_ave>p_3swa then deltap_3swa:=p_3_ave-p_3swa;
if p_1<p_1wa then deltap_1wa:=p_1wa-p_1;
else if p_1=p_1wa then deltap_1wa:=0;
else if p_1>p_1wa then deltap_1wa:=p_1-p_1wa;
if p_2<p_2wa then deltap_2wa:=p_2wa-p_2;
else if p_2=p_2wa then deltap_2wa:=0;
else if p_2>p_2wa then deltap_2wa:=p_2-p_2wa;
if p_4<p_4wa then deltap_4wa:=p_4wa-p_4;
else if p_4=p_4wa then deltap_4wa:=0;
else if p_4>p_4wa then deltap_4wa:=p_4-p_4wa;
if p_5<p_5lwa then deltap_5lwa:=p_5lwa-p_5;
else if p_5=p_5lwa then deltap_5lwa:=0;
else if p_5>p_5lwa then deltap_5lwa:=p_5-p_5lwa;
if p_5<p_5rwa then deltap_5rwa:=p_5rwa-p_5;
else if p_5=p_5rwa then deltap_5rwa:=0;
else if p_5>p_5rwa then deltap_5rwa:=p_5-p_5rwa;
if p_3_ave<p_1 then deltap_31:=p_1-p_3_ave;
else if p_3_ave=p_1 then deltap_31:=0;
else if p_3_ave>p_1 then deltap_31:=p_3_ave-p_1;
if p_3_ave<p_2 then deltap_32:=p_2-p_3_ave;
else if p_3_ave=p_2 then deltap_32:=0;
else if p_3_ave>p_2 then deltap_32:=p_3_ave-p_2;
if p_3_ave<p_4 then deltap_34:=p_4-p_3_ave;
else if p_3_ave=p_4 then deltap_34:=0;
else if p_3_ave>p_4 then deltap_34:=p_3_ave-p_4;
if p_3_ave<p_5 then deltap_35:=p_5-p_3_ave;
else if p_3_ave=p_5 then deltap_35:=0;
else if p_3_ave>p_5 then deltap_35:=p_3_ave-p_5;
if p_3_ave<p_6 then deltap_36:=p_6-p_3_ave;
else if p_3_ave=p_6 then deltap_36:=0;
else if p_3_ave>p_6 then deltap_36:=p_3_ave-p_6;
v_1a=c*abs(deltap_1wa)^n*d_0;
v_2a=c*abs(deltap_2wa)^n*d_0;
v_3swa=c*abs(deltap_3swa)^n*d_0;
v_3nda=c*abs(deltap_3nda)^n*d_1;
v_3nwa=c*abs(deltap_3nwa)^n*d_0;
v_4a=c*abs(deltap_4wa)^n*d_0;
v_5la=c*abs(deltap_5lwa)^n*d_0;
v_5ra=c*abs(deltap_5rwa)^n*d_0;
v_6a=c*abs(deltap_6a)^n*d_1;
v_31=c*abs(deltap_31)^n*d_2;
v_32=c*abs(deltap_32)^n*d_2;
v_34=c*abs(deltap_34)^n*d_2;
v_35=c*abs(deltap_35)^n*d_2;
v_36=c*abs(deltap_36)^n*d_2;
if p_1<p_1wa then v_1a_star:=v_1a;
else if p_1=p_1wa then v_1a_star:=0;
else if p_1>p_1wa then v_1a_star:=-v_1a;
if p_2<p_2wa then v_2a_star:=v_2a;
else if p_2=p_2wa then v_2a_star:=0;
else if p_2>p_2wa then v_2a_star:=-v_2a;
if p_3_ave<p_3nda then v_3nda_star:=v_3nda;
else if p_3_ave=p_3nda then v_3nda_star:=0;
else if p_3_ave>p_3nda then v_3nda_star:=-v_3nda;
if p_3_ave<p_3nwa then v_3nwa_star:=v_3nwa;
else if p_3_ave=p_3nwa then v_3nwa_star:=0;
else if p_3_ave>p_3nwa then v_3nwa_star:=-v_3nwa;
if p_3_ave<p_3swa then v_3swa_star:=v_3swa;
else if p_3_ave=p_3swa then v_3swa_star:=0;
else if p_3_ave>p_3swa then v_3swa_star:=-v_3swa;
if p_4<p_4wa then v_4a_star:=v_4a;
else if p_4=p_4wa then v_4a_star:=0;
else if p_4>p_4wa then v_4a_star:=-v_4a;
if p_5<p_5lwa then v_5la_star:=v_5la;
else if p_5=p_5lwa then v_5la_star:=0;
else if p_5>p_5lwa then v_5la_star:=-v_5la;
if p_5<p_5rwa then v_5ra_star:=v_5ra;
else if p_5=p_5rwa then v_5ra_star:=0;
else if p_5>p_5rwa then v_5ra_star:=-v_5ra;
if p_6<p_6sda then v_6a_star:=v_6a;
else if p_6=p_6sda then v_6a_star:=0;
else if p_6>p_6sda then v_6a_star:=-v_6a;
if p_1<p_3_ave then v_31_star:=v_31;
else if p_1=p_3_ave then v_31_star:=0;
else if p_1>p_3_ave then v_31_star:=-v_31;
if p_2<p_3_ave then v_32_star:=v_32;
else if p_2=p_3_ave then v_32_star:=0;
else if p_2>p_3_ave then v_32_star:=-v_32;
if p_4<p_3_ave then v_34_star:=v_34;
else if p_4=p_3_ave then v_34_star:=0;
else if p_4>p_3_ave then v_34_star:=-v_34;
if p_5<p_3_ave then v_35_star:=v_35;
else if p_5=p_3_ave then v_35_star:=0;
else if p_5>p_3_ave then v_35_star:=-v_35;
if p_6<p_3_ave then v_36_star:=v_36;
else if p_6=p_3_ave then v_36_star:=0;
else if p_6>p_3_ave then v_36_star:=-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;
SSRE_1a=sum(i=0:97)((deltap_1a_exp-deltap_1a)^2/deltap_1a_exp^2);
SSRE_2a=sum(i=0:97)((deltap_2a_exp-deltap_2a)^2/deltap_2a_exp^2);
SSRE_4a=sum(i=0:97)((deltap_4a_exp-deltap_4a)^2/deltap_4a_exp^2);
SSRE_5a=sum(i=0:97)((deltap_5a_exp-deltap_5a)^2/deltap_5a_exp^2);
SSRE_6a=sum(i=0:97)((deltap_6a_exp-deltap_6a)^2/deltap_6a_exp^2);
SSRE_3nda=sum(i=0:97)((deltap_3ndagauge_exp-deltap_3nda)^2/deltap_3ndagauge_exp^2);
SSRE_3nwa=sum(i=0:97)((deltap_3nwagauge_exp-deltap_3nwa)^2/deltap_3nwagauge_exp^2);
SSRE_3swa=sum(i=0:97)((deltap_3swagauge_exp-deltap_3swa)^2/deltap_3swagauge_exp^2);
SSRE_total=SSRE_1a+SSRE_2a+SSRE_4a+SSRE_5a+SSRE_6a+SSRE_3nda+SSRE_3nwa+SSRE_3swa;
minfunction SSRE_total;
页:
[1]