惩罚函数法

惩罚函数法

#include"stdio.h"

#include"stdlib.h"

#include"math.h"

const int kkg=3;

double r0;

double f(doublex[])

{doubleff;

ff=pow((x[0]-5),2)+4*pow((x[1]-6),2);

return(ff);

}

/*约束条件子程序*/

void strain(doublex[],doubleg[])

{g[0]=10-x[0];

g[1]=10+x[0]-x[1];

g[2]=25-pow((x[0]-8),2)+pow((x[1]-12),2);

}

/*惩罚函数子程序*/

double objf(doublep[])

{inti;

double ff,sg,*g;

g=(double*)malloc(kkg*sizeof(double));

sg=0;

strain(p,g);

for(i=0;i<kkg;i++)

{if(*(g+i)>0)

sg=sg+r0/(*(g+i));

else

sg=sg+r0*(1e+10);

}

free(g);

ff=f(p)+sg;

return(ff);

}

/*进退函数*/

void jtf(doublex0[],doubleh0,double s[],intn,double a[],doubleb[]){

int i;

double *xx[3],h,f1,f2,f3;

for (i=0;i<3;i++)

xx[i]=(double*)malloc(n*sizeof(double));

h=h0;

for(i=0;i<n;i++)

*(xx[0]+i)=x0[i];

f1=objf(xx[0]);

for(i=0;i<n;i++)

*(xx[1]+i)=*(xx[0]+i)+h*s[i];

f2=objf(xx[1]);

if(f2>=f1)

{h=-h0;

for(i=0;i<n;i++)

*(xx[2]+i)=*(xx[0]+i);

f3=f1;

for(i=0;i<n;i++)

{*(xx[0]+i)=*(xx[1]+i);

*(xx[1]+i)=*(xx[2]+i);

}

f1=f2;

f2=f3;

}

for(;;)

{h=2.0*h;

for(i=0;i<n;i++)

*(xx[2]+i)=*(xx[1]+i)+h*s[i];

f3=objf(xx[2]);

if(f2<f3)

break;

else

{for(i=0;i<n;i++)

{*(xx[0]+i)=*(xx[1]+i);

*(xx[1]+i)=*(xx[2]+i);

}

f1=f2;

f2=f3;

}

}

if(h<0)

for(i=0;i<n;i++)

{a[i]=*(xx[2]+i);

b[i]=*(xx[0]+i);

}

for(i=0;i<n;i++)

{a[i]=*(xx[0]+i);

b[i]=*(xx[2]+i);

}

for(i=0;i<3;i++)

free(xx[i]);

}

/*黄金分割程序*/

double gold (doublea[],doubleb[],doubleeps,int n,double x[])

{

int i;

double f1,f2,*xx[2],ff,q,w;

for(i=0;i<2;i++)

xx[i]=(double*)malloc(n*sizeof(double));

for(i=0;i<n;i++)

{*(xx[0]+i)=a[i]+0.618*(b[i]-a[i]);

*(xx[1]+i)=a[i]+0.382*(b[i]-a[i]);

}

f1=objf(xx[0]);

f2=objf(xx[1]);

{if(f1>f2)

{for(i=0;i<n;i++)

{b[i]=*(xx[0]+i);

*(xx[0]+i)=*(xx[1]+i);

}

f1=f2;

for(i=0;i<n;i++)

*(xx[1]+i)=a[i]+0.382*(b[i]-a[i]);

f2=objf(xx[1]);

}

else

{for(i=0;i<n;i++)

{a[i]=*(xx[1]+i);

*(xx[1]+i)=*(xx[0]+i);

}

f2=f1;

for(i=0;i<n;i++)

*(xx[0]+i)=a[i]+0.618*(b[i]-a[i]);

f1=objf(xx[0]);

}

q=0;

for(i=0;i<n;i++)

q=q+(b[i]-a[i])*(b[i]-a[i]);

w=sqrt(q);

}while(w>eps);

for(i=0;i<n;i++)

x[i]=0.5*(a[i]+b[i]);

ff=objf(x);

for(i=0;i<2;i++)

free(xx[i]);

return(ff);

}

double oneoptim(doublex0[],doubles[],doubleh0,double epsg,int n,double x[])

{doubleff,*a,*b;

a=(double*)malloc(n*sizeof(double));

b=(double*)malloc(n*sizeof(double));

jtf(x0,h0,s,n,a,b);

ff=gold(a,b,epsg,n,x);

free(a);

free(b);

return(ff);

}

/*powell子程序*/

double powell(doublep[],doubleh0,double eps,double epsg,int n,double x[]){inti,j,m;

double *xx[4],*s,*ss;

double f,f0,f1,f2,f3,fx,dlt,df,sdx,q,d;

ss=(double*)malloc(n*(n+1)*sizeof(double));

s=(double*)malloc(n*sizeof(doub

le));

for(i=0;i<=n;i++)

{for(j=0;j<=n;j++)

*(ss+i*(n+1)+j)=0;

*(ss+i*(n+1)+i)=1;

}

for(i=0;i<4;i++)

xx[i]=(double*)malloc(n*sizeof(double));for(i=0;i<n;i++)

*(xx[0]+i)=p[i];

for(;;)

{for(i=0;i<n;i++)

{*(xx[1]+i)=*(xx[0]+i);

x[i]=*(xx[1]+i);

}

f0=f1=objf(x);

dlt=-1;

for(j=0;j<n;j++)

{for(i=0;i<n;i++)

{*(xx[0]+i)=x[i];

*(s+i)=*(ss+i*(n+1)+j);

}

f=oneoptim(xx[0],s,h0,epsg,n,x);df=f0-f;

if(df>dlt)

{dlt=df;

m=j;

}

}

sdx=0;

for(i=0;i<n;i++)

sdx=sdx+fabs(x[i]-(*(xx[1]+i)));if(sdx<eps)

{

for(i=0;i<4;i++)

free(xx[i]);

return(f);

}

for(i=0;i<n;i++)

*(xx[2]+i)=x[i];

f2=f;

for(i=0;i<n;i++)

{*(xx[3]+i)=2*(*(xx[2]+i)-(*(xx[1]+i)));x[i]=*(xx[3]+i);

}

fx=objf(x);

f3=fx;

q=(f1-2*f2+f3)*(f1-f2-dlt)*(f1-f2-dlt);d=0.5*dlt*(f1-f3)*(f1-f3);

if((f3<f1)||(q<d))

{if(f2<=f3)

for(i=0;i<n;i++)

*(xx[0]+i)=*(xx[2]+i);

else

for(i=0;i<n;i++)

*(xx[0]+i)=*(xx[3]+i);

}

else

{for(i=0;i<n;i++)

{*(ss+(i+1)*(n+1))=x[i]-(*(xx[1]+i));

*(s+i)=*(ss+(i+1)*(n+1));

}

f=oneoptim(xx[0],s,h0,epsg,n,x);

for(i=0;i<n;i++)

*(xx[0]+i)=x[i];

for(j=m+1;j<=n;j++)

for(i=0;i<n;i++)

*(ss+i*(n+1)+j-1)=*(ss+i*(n+1)+j);

}

}

}

/*内点惩罚函数主程序*/

void main()

{int i;

double p[]={3,4};

double fom,fxo,c,x[2];

c=0.1;

r0=120;

fom=100;

do

{

fxo=powell(p,0.1,0.001,0.0001,2,x);

if(fabs(fom-fxo)>0.001)

{

fom=fxo;

r0=c*r0;

for(i=0;i<2;i++)

*(p+i)=x[i];

}

else

{

printf("输出最优点及其目标函数值：\n");

printf("x[0]=%f,x[1]=%f,ff=%f\n",x[0],x[1],fxo);return;

}

}while(1);

}

惩罚函数法

#include"stdio.h"

#include"stdlib.h"

#include"math.h"

const int kkg=3;

double r0;

double f(doublex[])

{doubleff;

ff=pow((x[0]-5),2)+4*pow((x[1]-6),2);

return(ff);

}

/*约束条件子程序*/

void strain(doublex[],doubleg[])

{g[0]=10-x[0];

g[1]=10+x[0]-x[1];

g[2]=25-pow((x[0]-8),2)+pow((x[1]-12),2);

}

/*惩罚函数子程序*/

double objf(doublep[])

{inti;

double ff,sg,*g;

g=(double*)malloc(kkg*sizeof(double));

sg=0;

strain(p,g);

for(i=0;i<kkg;i++)

{if(*(g+i)>0)

sg=sg+r0/(*(g+i));

else

sg=sg+r0*(1e+10);

}

free(g);

ff=f(p)+sg;

return(ff);

}

/*进退函数*/

void jtf(doublex0[],doubleh0,double s[],intn,double a[],doubleb[]){

int i;

double *xx[3],h,f1,f2,f3;

for (i=0;i<3;i++)

xx[i]=(double*)malloc(n*sizeof(double));

h=h0;

for(i=0;i<n;i++)

*(xx[0]+i)=x0[i];

f1=objf(xx[0]);

for(i=0;i<n;i++)

*(xx[1]+i)=*(xx[0]+i)+h*s[i];

f2=objf(xx[1]);

if(f2>=f1)

{h=-h0;

for(i=0;i<n;i++)

*(xx[2]+i)=*(xx[0]+i);

f3=f1;

for(i=0;i<n;i++)

{*(xx[0]+i)=*(xx[1]+i);

*(xx[1]+i)=*(xx[2]+i);

}

f1=f2;

f2=f3;

}

for(;;)

{h=2.0*h;

for(i=0;i<n;i++)

*(xx[2]+i)=*(xx[1]+i)+h*s[i];

f3=objf(xx[2]);

if(f2<f3)

break;

else

{for(i=0;i<n;i++)

{*(xx[0]+i)=*(xx[1]+i);

*(xx[1]+i)=*(xx[2]+i);

}

f1=f2;

f2=f3;

}

}

if(h<0)

for(i=0;i<n;i++)

{a[i]=*(xx[2]+i);

b[i]=*(xx[0]+i);

}

for(i=0;i<n;i++)

{a[i]=*(xx[0]+i);

b[i]=*(xx[2]+i);

}

for(i=0;i<3;i++)

free(xx[i]);

}

/*黄金分割程序*/

double gold (doublea[],doubleb[],doubleeps,int n,double x[])

{

int i;

double f1,f2,*xx[2],ff,q,w;

for(i=0;i<2;i++)

xx[i]=(double*)malloc(n*sizeof(double));

for(i=0;i<n;i++)

{*(xx[0]+i)=a[i]+0.618*(b[i]-a[i]);

*(xx[1]+i)=a[i]+0.382*(b[i]-a[i]);

}

f1=objf(xx[0]);

f2=objf(xx[1]);

{if(f1>f2)

{for(i=0;i<n;i++)

{b[i]=*(xx[0]+i);

*(xx[0]+i)=*(xx[1]+i);

}

f1=f2;

for(i=0;i<n;i++)

*(xx[1]+i)=a[i]+0.382*(b[i]-a[i]);

f2=objf(xx[1]);

}

else

{for(i=0;i<n;i++)

{a[i]=*(xx[1]+i);

*(xx[1]+i)=*(xx[0]+i);

}

f2=f1;

for(i=0;i<n;i++)

*(xx[0]+i)=a[i]+0.618*(b[i]-a[i]);

f1=objf(xx[0]);

}

q=0;

for(i=0;i<n;i++)

q=q+(b[i]-a[i])*(b[i]-a[i]);

w=sqrt(q);

}while(w>eps);

for(i=0;i<n;i++)

x[i]=0.5*(a[i]+b[i]);

ff=objf(x);

for(i=0;i<2;i++)

free(xx[i]);

return(ff);

}

double oneoptim(doublex0[],doubles[],doubleh0,double epsg,int n,double x[])

{doubleff,*a,*b;

a=(double*)malloc(n*sizeof(double));

b=(double*)malloc(n*sizeof(double));

jtf(x0,h0,s,n,a,b);

ff=gold(a,b,epsg,n,x);

free(a);

free(b);

return(ff);

}

/*powell子程序*/

double powell(doublep[],doubleh0,double eps,double epsg,int n,double x[]){inti,j,m;

double *xx[4],*s,*ss;

double f,f0,f1,f2,f3,fx,dlt,df,sdx,q,d;

ss=(double*)malloc(n*(n+1)*sizeof(double));

s=(double*)malloc(n*sizeof(doub

le));

for(i=0;i<=n;i++)

{for(j=0;j<=n;j++)

*(ss+i*(n+1)+j)=0;

*(ss+i*(n+1)+i)=1;

}

for(i=0;i<4;i++)

xx[i]=(double*)malloc(n*sizeof(double));for(i=0;i<n;i++)

*(xx[0]+i)=p[i];

for(;;)

{for(i=0;i<n;i++)

{*(xx[1]+i)=*(xx[0]+i);

x[i]=*(xx[1]+i);

}

f0=f1=objf(x);

dlt=-1;

for(j=0;j<n;j++)

{for(i=0;i<n;i++)

{*(xx[0]+i)=x[i];

*(s+i)=*(ss+i*(n+1)+j);

}

f=oneoptim(xx[0],s,h0,epsg,n,x);df=f0-f;

if(df>dlt)

{dlt=df;

m=j;

}

}

sdx=0;

for(i=0;i<n;i++)

sdx=sdx+fabs(x[i]-(*(xx[1]+i)));if(sdx<eps)

{

for(i=0;i<4;i++)

free(xx[i]);

return(f);

}

for(i=0;i<n;i++)

*(xx[2]+i)=x[i];

f2=f;

for(i=0;i<n;i++)

{*(xx[3]+i)=2*(*(xx[2]+i)-(*(xx[1]+i)));x[i]=*(xx[3]+i);

}

fx=objf(x);

f3=fx;

q=(f1-2*f2+f3)*(f1-f2-dlt)*(f1-f2-dlt);d=0.5*dlt*(f1-f3)*(f1-f3);

if((f3<f1)||(q<d))

{if(f2<=f3)

for(i=0;i<n;i++)

*(xx[0]+i)=*(xx[2]+i);

else

for(i=0;i<n;i++)

*(xx[0]+i)=*(xx[3]+i);

}

else

{for(i=0;i<n;i++)

{*(ss+(i+1)*(n+1))=x[i]-(*(xx[1]+i));

*(s+i)=*(ss+(i+1)*(n+1));

}

f=oneoptim(xx[0],s,h0,epsg,n,x);

for(i=0;i<n;i++)

*(xx[0]+i)=x[i];

for(j=m+1;j<=n;j++)

for(i=0;i<n;i++)

*(ss+i*(n+1)+j-1)=*(ss+i*(n+1)+j);

}

}

}

/*内点惩罚函数主程序*/

void main()

{int i;

double p[]={3,4};

double fom,fxo,c,x[2];

c=0.1;

r0=120;

fom=100;

do

{

fxo=powell(p,0.1,0.001,0.0001,2,x);

if(fabs(fom-fxo)>0.001)

{

fom=fxo;

r0=c*r0;

for(i=0;i<2;i++)

*(p+i)=x[i];

}

else

{

printf("输出最优点及其目标函数值：\n");

printf("x[0]=%f,x[1]=%f,ff=%f\n",x[0],x[1],fxo);return;

}

}while(1);

}


相关文章

  • 基于VB外点惩罚函数法的实现

    • Equipment Manufactring Technology No.3,2011 基于VB外点惩罚函数法的实现 殷晓飞1,任晓丹2 (1．呼和浩特职业学院机电工程学院,内蒙古呼和浩特010051: 2．内蒙古机电职业技术学院电气工程系 ...查看


  • 求解非线性规划问题的遗传算法设计与实现

    • 摘 要 非线性规划在工程.管理.经济.科研.军事等方面都有广泛的应用.传统的解决非线性规划问题的方法,如梯度法.罚函数法.拉格朗日乘子法等,稳定性差,对函数初值和函数性态要求较高,且容易陷入局部最优解. 遗传算法是模拟达尔文的遗传选择和自然 ...查看


  • 型钢混凝土矩形截面梁优化设计

    • 第39卷增刊2 20 哈尔滨工业大学学报 JOURNALOFHARBININSTITUTEOF VoL39 Sup・2 07年8月 TECHNOLOGY Aug・2007 型钢混凝土矩形截面梁优化设计 曾 磊,郑山锁,查春光,李磊,邓国专, ...查看


  • 运筹学外点法

    • (一) 实验目的 熟练掌握外点法.内点法原理并可以在matlab 熟练运行. (二) 问题描述 目标函数:z =min x 1+x 2 s.t. –x 21+x 2≥0 x 1≥0 (三) 算法介绍 外点法: 对于混合约束问题 min f( ...查看


  • 优化设计 孙靖民 课后答案第6章习题解答

    • 第六章习题解答 1． 已知约束优化问题: min f (x ) =(x 1-2) 2+(x 2-1) 2 s ⋅t g 1(x ) =x 12-x 2≤0g 2(x ) =x 1+x 2-2≤0 试从第k 次的迭代点x (k ) =[-12 ...查看


  • 惩罚函数的内点法

    • 201-2014 (1) 专业课程实践论文 内点法 一.算法理论 内点法总是从可行域的内点出发,并保持在可行域内进行搜索,因此这种方法适用于只有不等式约束条件的问题 内点法据图计算步骤: 1．给定初x(0)∈intD,允许误差ε〉0,初始参 ...查看


  • 修正免疫克隆约束多目标优化算法

    • 软件学报ISSN 1000-9825, CODEN RUXUEW E-mail: [email protected] Journal of Software,2012,23(7):1773−1786 [doi: 10.3724/SP.J.100 ...查看


  • 基于改进遗传算法的汽车混流装配线物料配送路径优化

    • 物流科技2016年第2期 k矛stics Sci-Tech No．2,2016 ・理论研究・ 文章编号:1002-3100(2016)02-0016-05 基于改进遗传算法的汔车混流装配线物料配送路径优化 0删．吼iza60n ofMate ...查看


  • 马尔可夫决策过程实例讲解

    • Machine Learning 16-Reinforcement Learning 之前我们学过3个部分的内容:监督学习.学习理论.半监督学习.现在我们来学习第四部分:自增强学习. 在监督学习中,给定了训练集以及对应的标签y ,算法要做的 ...查看


热门内容