蒙特卡洛方法及数学建模算法与应用第二章习题

蒙特卡洛的方法又称为随机取样法，所以他是根据大量随机数来跑出满意解

比如上一章用遗传算法的完全可以用蒙特卡洛方法（一口血吐出....）

 以下为c语言代码：
#include <stdio.h>
#include <stdlib.h>
#include <time.h>
#include <string.h>
int check(int x[])
{
if (x[0]+x[1]+x[2]+x[3]+x[4]>400)
return 0;
if (x[0]+2*x[1]+2*x[2]+x[3]+6*x[4]>800)
return 0;
if (2*x[0]+x[1]+6*x[2]>200)
return 0;
if (x[2]+x[3]+5*x[4]>200)
return 0;
return 1;

}
int main()
{
int x[5];
int ans[5];
int i,j,max;
srand(time(NULL));
max=-2;
memset(ans,0,sizeof(ans));
for (i=0;i<100000;i++)
{
for (j=0;j<5;j++)
x[j]=rand()%100;
if (check(x))
if ((x[0]*x[0]+x[1]*x[1]+x[2]*x[2]*3+4*x[3]*x[3]+2*x[4]*x[4]-8*x[0]-2*x[1]-3*x[2]-x[3]-2*x[4])>max)
{
for (j=0;j<5;j++)
ans[j]=x[j];
max=x[0]*x[0]+x[1]*x[1]+x[2]*x[2]*3+4*x[3]*x[3]+2*x[4]*x[4]-8*x[0]-2*x[1]-3*x[2]-x[3]-2*x[4];
}
}
printf("最大值为:%d\n",max);
for (j=0;j<5;j++)
printf("%d ",ans[j]);
}


2.1：
换元，另y=x1*x2，接下来挺难想到的~~

 0<=x1<=1;

0<=x2<=1;

故(x1-1)*(x2-1)>=0

   x1*x2-x1-x2+1>=0

   x1*x2>=x1+x2-1

   y>=x1+x2-1

   (x1-1)*x2<=0

  x1*x2-x2<=0

   y<=x2

   同理得y<=x1

  综上       x1+x2-1<=y<=x1

                 x1+x2-1<=y<=x2

     >> f=[-1 0 1 -1];
>> a=[-2 3 1 0;
      1 1 0 -1;
      -1 0 0 1;
      0 -1 0 1];
>> b=[3;1;0;0];
>> c=[1,2,3,4];
>> [x,fval]=intlinprog(f,c,a,b,[],[],zeros(4,1),ones(4,1))


2.2
   设xi为第i个小区是否被选，选为1,不被选为0

  minz=x1+x2+x3+x4+x5+x6

   对于每一个小区又必须要有一所学校

 x1+x2+x3>=1

x2+x4>=1

x3+x5>=1

x4+x6>=1

x1+x2+x3>=1

x5+x6>=1

x1>=1

x2+x4+x6>=1

>> f=[1,1,1,1,1,1];
>> a=[1,2,3,4,5,6];
>> b=[-1 -1 -1 0 0 0;0 -1 0 -1 0 0;0 0 -1 0 -1 0;0 0 0 -1 0 -1;-1 -1 -1 0 0 0;0 0 0 0 -1 -1;-1 0 0 0 0 0;0 -1 0 -1 0 -1];
>> c=[-1;-1;-1;-1;-1;-1;-1;-1];
>> [x,fval]=intlinprog(f,a,b,c,[],[],zeros(6,1),ones(6,1))

例2.3题目描述貌似有点问题，搜了下答案应该意思是每个企业有且只有1台，每一个设备至少要购入一件。
整数规划问题，用蒙特卡洛方法可解。

#include <stdio.h>
#include <stdlib.h>
#include <time.h>
#include <string.h>
int value[7][5];
void init()
{
value[1][1]=4;
value[1][2]=2;
value[1][3]=3;
value[1][4]=4;
value[2][1]=6;
value[2][2]=4;
value[2][3]=5;
value[2][4]=5;
value[3][1]=7;
value[3][2]=6;
value[3][3]=7;
value[3][4]=6;
value[4][1]=7;
value[4][2]=8;
value[4][3]=8;
value[4][4]=6;
value[5][1]=7;
value[5][2]=9;
value[5][3]=8;
value[5][4]=6;
value[6][1]=7;
value[6][2]=10;
value[6][3]=8;
value[6][4]=6;
}
int check(int x[])
{
int i,j;
int vis[5];
memset(vis,0,sizeof(vis));
for (i=1;i<=6;i++)
vis[x[i]]++;
for (i=1;i<=4;i++)
if (vis[i]==0)
return 0;
return 1;
}
int find(int x[])
{
int sum,i;
sum=0;
for (i=1;i<=6;i++)
sum=sum+value[i][x[i]];
return sum;
}
int main()
{
int i,j,k;
int x[7];
int ans[7];
int max;
srand(time(NULL));
init();
max=-1;
memset(ans,0,sizeof(ans));
for (i=1;i<=10000;i++)
{
for (j=1;j<=6;j++)
x[j]=rand()%4+1;
if (check(x))
if (find(x)>max)
{
max=find(x);
for (j=1;j<=6;j++)
ans[j]=x[j];
}
}
printf("最大值为:%d\n",max);
for (i=1;i<=6;i++)
printf("%d ",ans[i]);

}
最大利润为44