分水系统模型.cpp

解线性规划

          break;   }                               //break必须有,达到前两条件后调处,k不再自加
	}

	void main()
	{   
   int i=0,j=0,k,H=0,I=1,J=0,II;
   float C[n][n],A[n],B[n],h1,h2,h3,FDC=10;
   float mu=100;
   float X0[n][n],X[n][n];
   float a[n],b[n];

   double f,f0;
	float   t,sum=0,p[n][n],p0[n][n],g[n][n],g0[n][n],rfa,e=1e-3,e1,e2,e3;
   

   FILE *fp1,*fp2,*fp3;

   for(i=0;i<n;i++)
	   for(j=0;j<n;j++)
		   X0[i][j]=1;         // X 被初始化 均为1
  fp1=fopen("参数库.txt","w+");
  fprintf(fp1,"2.5 1.60 4.55 2.18 2.18 10000 589\n");                        //题目参量矩阵
  fprintf(fp1,"2.08 0 2.56 1.40 0.25 10000 49\n");
  fprintf(fp1,"0.95 0 1.53 0.15 0.12 10000 295\n");
  fprintf(fp1,"2.46 0 2.94 1.66 1.28 10000 48\n");
  fprintf(fp1,"1.11 0.099 1.87 0.42 0.42 10000 5.4\n");
  fprintf(fp1,"0.80 0 0.74 0 0 10000 18\n");
  
  //fprintf(fp1,"10000 10000 10000 10000 10000 10000 10000 10000 10000 10000 0\n");
  //fprintf(fp1,"10000 10000 10000 10000 10000 10000 10000 10000 10000 10000 0\n");
  //fprintf(fp1,"10000 10000 10000 10000 10000 10000 10000 10000 10000 10000 0\n");
  //fprintf(fp1,"10000 10000 10000 10000 10000 10000 10000 10000 10000 10000 0\n");
  fprintf(fp1,"486 97 90 295 36.4 0\n");//第11行:B[j]
  
   //fprintf(fp1,"1 2 3 4 5 6 7 8 9 0 10\n");                                        //测试矩阵
   //fprintf(fp1,"11 12 13 14 15 16 17 18 19 0 10\n");
   //fprintf(fp1,"21 22 23 24 25 26 27 28 29 0 10\n");
   //fprintf(fp1,"31 32 33 34 35 36 37 38 39 0 10\n");
   //fprintf(fp1,"41 42 43 44 45 46 47 48 49 0 10\n");
   //fprintf(fp1,"51 52 53 54 55 56 57 58 59 0 10\n");
   //fprintf(fp1,"61 62 63 64 65 66 67 68 69 0 10\n");
   //fprintf(fp1,"71 72 73 74 75 76 77 78 79 0 10\n");
   //fprintf(fp1,"81 82 83 84 85 86 87 88 89 0 10\n");
   //fprintf(fp1,"91 92 93 94 95 96 97 98 99 0 10\n");
   //fprintf(fp1,"10 10 10 10 10 10 10 10 10 10\n");//第11行:B[j]
   
  fclose(fp1);

   fp1=fopen("参数库.txt","r");
   for(i=0;i<n;i++)
   { for(j=0;j<n;j++)     fscanf(fp1,"%f",&C[i][j]);//输入C[i][j]
   fscanf(fp1,"%f",&A[i]); }                      //输入A[i][j]——行和
   for(i=0;i<n;i++)   fscanf(fp1,"%f",&B[i]);       //输入B[i][j]——列和
   fclose(fp1);
   
   fp2=fopen("结果库.txt","w+");//输出文件,格式(规整)化系数矩阵
   fprintf(fp2,"A====================================A\n");
   for(i=0;i<n;i++)
   {   for(j=0;j<n;j++)
      fprintf(fp2,"%5.3f\t",C[i][j]);
      fprintf(fp2,"%-5.3f\t",A[i]);
      fprintf(fp2,"\n"); }
   fprintf(fp2,"\n");
      for(i=0;i<n;i++)   
      fprintf(fp2,"%8.3f",B[i]);
   fprintf(fp2,"\nA===========A===========A============A\n");
   fclose(fp2);
   

for(II=0;II<20;II++)                                           //最外层_起始 

{









   k=1;  J=0;  j=0;                                         //进入外罚法======================================

   fp3=fopen("ls_共厄梯度解目标.txt","w+");
       fprintf(fp3,"START=========================================================================WRITE\n",I++);
   fclose(fp3); ++I;
ceshihanghe(a,X0);
ceshiliehe(a,X0);
fp3=fopen("ls_共厄梯度解目标.txt","a");
   fprintf(fp3,"\t=====mu=%-17.e===========\n",mu);               
	   fclose(fp3);
   eq(X,X0);


   f=F_min(C,X0,A,B,mu);
   F_DF(g,C,X0,A,B,mu) ;

   
   fp3=fopen("ls_共厄梯度解目标.txt","a");
       fprintf(fp3,"\nFx=%-17.5e\n",F_min(C,X,A,B,mu));
       rf(C,X,A,B);

   fclose(fp3);
   
   
   do
   {    k=0;     uni_eq(p0,g);

   
   
	for(;H!=1&&k<300;k++)
	{  eq(X0,X);
	   f0=f;
	   eq(g0,g);

rf(C,X,A,B);


	   ls(p0,X0,C,X,A,B,mu);     
	   //x0 has been put into x
		                     
	   fp3=fopen("ls_共厄梯度解目标.txt","a");
       fprintf(fp3,"ls%d=====================================================================ls%d\n",J,J);
       rf(C,X,A,B);
	   for(i=0;i<n;i++)
	   { for(j=0;j<n;j++)
       fprintf(fp3,"x%-d%d=%-17.5e",i,j,X[i][j]);
	   fprintf(fp3,"\n"); }
       fclose(fp3); 
	   ++J;
	   f=F_min(C,X,A,B,mu);
	   F_DF(g,C,X,A,B,mu) ;
ceshihanghe(a,X0);
ceshiliehe(a,X0);

	   e1=e2=1e-5;e3=1e-4;                    //H order
	   h1=h2=h3=0;
	      for(i=0;i<n;i++)
			  for(j=0;j<n;j++)
			  { h1=h1+fabs(g[i][j]);
	          	h2=h2+fabs(X[i][j]-X0[i][j]);
	        	h3=h3+fabs(X0[i][j]);  }
		h2=h2/(h3+1);
		double(h3)=fabs(F_min(C,X,A,B,mu)-F_min(C,X0,A,B,mu))/(fabs(F_min(C,X0,A,B,mu))+1);
	    fp3=fopen("ls_共厄梯度解目标.txt","a");
		fprintf(fp3,"\n\nh1=%-17.eh2=%-17.eh3=%-17.e",h1,h2,h3);
		fprintf(fp3,"\ne3=%-17.ee1=%-17.ee2=%-17.e\n",e3,e1,e2);
	     fclose(fp3);
		if(h1<e3&&h2<=e1&&h3<=e2)
		{ H=1;
	     fp3=fopen("ls_共厄梯度解目标.txt","a");
	     fprintf(fp3,"\n=======================H order has been passed.=======================\n");
	     for(i=0;i<n;i++)
			 for(j=0;j<n;j++)    
				 fprintf(fp3,"x%-d%-d=%-17.5e",i,j,X[i][j]);
	     fprintf(fp3,"\nH=%-18dFx=%-17.5ek=%-18dJ=%-17d\n",H,F_min(C,X,A,B,mu),k,J);
	     fprintf(fp3,"\nh1=%-17.eh2=%-17.eh3=%-17.10e\n",h1,h2,h3);
	     fprintf(fp3,"\ne3=%-17.ee1=%-17.ee2=%-17.e\n",e3,e1,e2);
	     fclose(fp3);
		break;   }
		fp3=fopen("ls_共厄梯度解目标.txt","a");
		fprintf(fp3,"\n");
	     fprintf(fp3,"H=%-18dFx=%-17.5ek=%-18dJ=%-17d\n",H,F_min(C,X,A,B,mu),k,J);
	     fprintf(fp3,"\t\tJump over the H order and go down.\n");
	     fclose(fp3);
	   if(k==n*n) { fp3=fopen("ls_共厄梯度解目标.txt","a");
	              fprintf(fp3,"\n======================从k=n跳出====================");
				  fclose(fp3);
				  break;}         //break from 'for' and continue 'do while'
	      rfa=product(g,g)/product(g0,g0);
	      add(-1,g,rfa,p0,p);               //p=-g+a*p0
	   if(fabs(product(p,g))<=e)  {fp3=fopen("ls_共厄梯度解目标.txt","a");
	                               fprintf(fp3,"\n======================从||p*g||<=e跳出====================");
				                   fclose(fp3);
								   break; }
	   else if(product(p,g)<-1*e)
	      eq(p0,p);
	      else  uni_eq(p0,p);
	     fp3=fopen("ls_共厄梯度解目标.txt","a");
	                               fprintf(fp3,"\n======================从k=k+1回for循环====================\n\n");               
								   fclose(fp3); }  //end for

   J++;
   fp3=fopen("ls_共厄梯度解目标.txt","a");
	  fprintf(fp3,"\n======================从while返回,继续while循环====================\n\n");
         fclose(fp3); 

   }while(H==0&&J<500);

                               //output on the square-eye
   for(i=0;i<n;i++)
   {		   for(j=0;j<n;j++)
   printf("x%-d%-d=%-17.5e",i,j,X[i][j]);
   printf("\nFx=%.5f\tk=%d\tj=%d\n",F_min(C,X,A,B,mu),k,j);         }



   if(log10(mu)+log10(rf(C,X,A,B))<-4) {fp3=fopen("ls_共厄梯度解目标.txt","a");
	                               fprintf(fp3,"\t=============计算终止=============mu=%-17.e======计算终止==========II=%d==============\n",mu,II);               
								   fclose(fp3);  break;}

   else { fp3=fopen("ls_共厄梯度解目标.txt","a");
	                               fprintf(fp3,"\t===计算没有终止==mu=%-17.e===========\n",mu);               
								   fclose(fp3);
								   mu=mu*10; eq(X0,X); printf("mu=%17.e\n",mu);J=0;}  


   }                                      //最外层_终止
   //改F_MIN为double

   

}