内点法

一、目的
使用目前熟悉的C 语言进行编程，实现所学的内点法对约束非线性规划问题进行求解，并正确运算测试结果。

二、内点惩罚函数法
内点法是从可行域内某一初始内点出发，在可行域内进行迭代的序列极小化方法。

它仅用于求解不等式约束优化问题。

三、问题的提出
用c 语言实现内点法的编程，并求解下面的非线性规划问题。

12
22121min ()..0
f x x x s t x x x =+-+≤≥
四、算法基础
内点法具体计算步骤为：
1、 给定初始点(0)i n t x D ∈，允许误差0ε>，初始参数10r >，缩小系数(0,1),1k β∈=；
2、 以(1)k x -为初始点，求解问题
min ()()
..int k f x r B x s t x D +∈
设求得极小点为()k x ；
3、 若()()k k r B x ε<则停，得近似解()k x ；否则令1,1k k r r k k β+==+，回2。

五、具体实现程序
#include <stdio.h>
#include <math.h>
#define MAXDEM 10
//This is the optional model:
struct MODEL
{
double (*g[10])(double *); //The object function
double (*f)(double *); //The subject functions,the max number is 10
int stNum; //The number of subject functions
} Model;
//******A example to get the optional value of function x1+x2
double f(double x[])
{
return x[0]*x[0]+x[1]*x[1];
}
//****下面的三个函数是不等式约束函数，要求约束函数>0
double g1(double x[])
{
return 0; //-x1*x1+x2>0
}
double g2(double x[])
{
return x[0]-1; //x1>0
}
//*****The function created for this method
double fun(double x[],double r)
{
int i;
double (*g)(double *);
double ret=0;
for(i=0;i<Model.stNum;i++)
{
g=Model.g[i];
if(g(x)<1e-10)
ret+=1e+30;
else
ret+=r/g(x);
}
return f(x)+ret;
}
//**********The following two functions are to get the optional value of fun
// in dimensions_turn method
void Search(double(*fun)(double *, double ),double Accur,int s,double X[],double r,int Dim) {
double S[MAXDEM],AidX[MAXDEM],h=1;
int i;
for(i=0;i<Dim;i++)
{
if(i==s)
S[i]=1;
else
S[i]=0;
}
while(fabs(h)>Accur)
{
for(i=0;i<Dim;i++)
AidX[i]=X[i]-h*S[i];
if(fun(X,r)>fun(AidX,r))
{
for(i=0;i<Dim;i++)
X[i]=AidX[i];
h*=2;
}
else
h=h*(-1.0/4.0);
}
}
void Demensions_Turn(double(*fun)(double *, double),double Accur,double* X,double r,int Dim) {
double AidX[MAXDEM],Mod=1;
int i;
while(Mod>Accur)
{
Mod=0;
for(i=0;i<Dim;i++)
AidX[i]=X[i];
for(i=0;i<Dim;i++)
Search(fun,Accur,i,X,r,Dim);
for(i=0;i<Dim;i++)
Mod+=(AidX[i]-X[i])*(AidX[i]-X[i]);
Mod=sqrt(Mod);
}
}
//The function is to get the optional value of the model
void Sum_In(double x[],double Accur,int Dim)
{
double Aidx[MAXDEM],r=1,mod=10;
int i;
printf("\tr");
for(i=0;i<Dim;i++)
printf("\tx%d",i+1);
printf("\tf(x*)\n");
while(mod>Accur&&r>1e-10)
{
mod=0;
r *=0.02;
for(i=0;i<Dim;i++)
Aidx[i]=x[i];
Demensions_Turn(fun,Accur,x,r,Dim);
printf("%10.9lf",r);
for(i=0;i<Dim;i++)
{
printf("%10.5lf",x[i]);
if(fabs(Aidx[i]-x[i])>mod)
mod=fabs(Aidx[i]-x[i]);
}
printf("%10.5lf",f(x));
printf("\n");
}
printf("\tr");
printf("\tx%d",i-1);
printf("\tx%d",i);
printf("\tf(x*)\n");
printf("%10.9lf",r);
printf("%10.5lf",x[i-2]);
printf("%10.5lf",x[i-1]);
printf("%10.5lf",f(x));
printf("\n");
}
int main()
{
double x[2]={1,0},Accur=0.000000005;
//Initial the Model
Model.f=f;
Model.g[0]=g1;
Model.g[1]=g2;
Model.stNum=2;
Sum_In(x,Accur,2);
return 0;
}
六、实验结果及分析
从运行结果可以看出，最优解为：*(0,0)T
x=
最优值为：
()0
f x=。