变尺度法的C语言实现

1、 变尺度法的C+语言实现 -基于成都理工大学最优化教材#include <iostream.h>#include "Matrix.h"#include<LIMITS.H>#define MAX_M 2048#include<math.h>#define beta (sqrt(5)-1)/2#define eclipse 0.19CMatrix hesse()CMatrix temp(2,2); temp.A00 = 8;temp.A01 = 0;temp.A10 = 0;temp.A11 = 2;return temp;double fu

2、n(double t,CMatrix &X1,CMatrix &X2)double x1 = X1.A00 + t*X2.A00;double x2 = X1.A10 + t*X2.A10;return 4*(x1-5)*(x1-5) + (x2-6)*(x2-6);double Max_Num(double a,double b)if(a>b)return a;elsereturn b;double Min_Num(double a,double b)if(a<b)return a;elsereturn b;void StepAdding_Search(doubl

3、e &a,double &b,CMatrix &mt1,CMatrix &mt2)double tMAX_M=0;double hMAX_M=0;double fMAX_M=0;double result=0;double p=2.1;t0=0;h0=1;f0=fun(t0,mt1,mt2);for(int k=0;k<MAX_M-1;k+)tk+1=tk+hk;fk+1=fun(tk+1,mt1,mt2);if(fk+1<fk)hk+1=p*hk;result=tk;tk=tk+1;fk=fk+1;else if(k = 0)hk=-hk;resu

4、lt=tk+1;elsea=Min_Num(result,tk+1);b=Max_Num(result,tk+1);double Fibonacci(double &a,double &b,CMatrix mt1,CMatrix mt2)double t2,t1,result_1,result_2;t2 = a + beta*(b-a);result_2 = fun(t2,mt1,mt1);while(true)t1 = a+b-t2;result_1 = fun(t1,mt1,mt2);if(fabs(t1-t2)<eclipse)return (t1+t2)/2;el

5、se if(result_1 < result_2)b = t1;t1 = t2;result_1 = result_2;t2 = a+beta*(b-a);result_2 = fun(t2,mt1,mt2);else if(result_1 = result_2)a = t1;b = t2;result_2 = fun(a,mt1,mt2);elsea = t2;t2 = t1;result_2 = result_1;double distance(CMatrix &X1,CMatrix &X2)double x1 = X1.A00 - X2.A00;double x

6、2 = X1.A10 - X2.A10;return x1*x1 + x2*x2;CMatrix diff_fun(CMatrix &mt)CMatrix temp(2,1);temp.A00 = 8*(mt.A00 - 5);temp.A10 = 2*(mt.A10 - 6);return temp;double fanshu(CMatrix mt)return sqrt(mt.A00*mt.A00 + mt.A10*mt.A10);void main()int i =1;int n = 2;double a ,b ;CMatrix temp_X1;CMatrix X1(2,1);C

7、Matrix H0(2,2);H0.A00 = 1;H0.A10 = 0;H0.A01 = 0;H0.A11 = 1;X1.A00 = 8;X1.A10 = 9;CMatrix m_hesse = hesse();CMatrix m_temp1 = (diff_fun(X1).Transpose()*diff_fun(X1);CMatrix m_temp2 = (diff_fun(X1).Transpose()*(m_hesse*diff_fun(X1);CMatrix X2 = X1 - diff_fun(X1).Matrix_Multiply(m_temp1.A00/m_temp2.A00

8、);CMatrix S0 = X2 - X1;CMatrix Y0 = diff_fun(X2) - diff_fun(X1);double m1,m2;m1 = (S0.Transpose() * Y0).A00;m2 = (Y0.Transpose()*H0*Y0).A00;CMatrix H1 = H0 + (S0*S0.Transpose().Matrix_Multiply(1/m1) + (H0*Y0*Y0.Transpose()*H0).Matrix_Multiply(1/m2);CMatrix m_search = (H1*diff_fun(X2).Matrix_Multiply

9、(-1);StepAdding_Search(a,b,X1,m_search);double m = Fibonacci(a,b,X1,m_search);CMatrix temp_X2 = X2;CMatrix X3 = X2 + m_search.Matrix_Multiply(m);CMatrix temp_H = H1;while(distance(temp_X2,X3) >eclipse)X2 = X3;if(i%2 = 1)m_temp1 = (diff_fun(X2).Transpose()*diff_fun(X2);m_temp2 = (diff_fun(X2).Tran

10、spose()*(m_hesse*diff_fun(X2);temp_X2 = X3;X3 = temp_X2 - diff_fun(X2).Matrix_Multiply(m_temp1.A00/m_temp2.A00);elseS0 = X3 - temp_X2;Y0 = diff_fun(X3) - diff_fun(temp_X2);m1 = (S0.Transpose() * Y0).A00;H0 = temp_H;m2 = (Y0.Transpose()*H0*Y0).A00;CMatrix H2 = temp_H + (S0*S0.Transpose().Matrix_Multiply(1/m1) + (H0*Y0*Y0.Transpose()*H0).Matrix_Multiply(1/m2);temp_H = H2;m_search = (H2*diff_fun(X2).Matrix_Multiply(-1);m = Fibonacci(a,b,X1,m_search);temp_X2 = X3;X3 = temp_X2 + m_search.Matrix_Multiply(m);/X2.print_Matrix();i+;cout<<"求解的结果是:"<<endl;cout<<"(&