蒙特卡洛方法解非线性规划问题

蒙特卡洛算法定义: 当所求解问题是某种随机事件出现的概率，或者是某个随机变量的期望值时，通过某种“实验”的方法，以这种事件出现的频率估计这一随机事件的概率，或者得到这个随机变量的某些数字特征，并将其作为问题的解。详情蒙特卡洛算法-百度百科。

这里我们将用蒙特卡洛方法求解一个非线性规划问题（以下简称NLP），NLP问题没有通用的解法，用Lingo解决这类问题是可行的，但是当运算量增大时，Lingo运行的效率并不高，往往花费大量时间在迭代上。有关NLP问题更多请见非线性规划-百度百科。

让我们想像一下这个问题，一个平面上有三个圆，那么如果用一条绳子将它们连接起来，那么这条绳子最短长度是多少呢？

首先我们给出三个圆在直角坐标系下的方程如下：

然后，哦！说好的绳子呢？别着急，这根绳子是这样一个约束，即一个过三个圆的三角形的周长，我们的目的是使这个三角形的周长最短。嗯！那么怎样实现呢？我想直角坐标系下两点的距离公式大家应该都记得，没错，就是它！因此，我们得到的目标函数是这个样子的：

是不是很简单？那么让我们求一下这个最小值吧，当然，最开始我将先用Lingo演示一下这个NLP问题的求法，程序如下：

min=@sqrt((x1-x2)^2+(y1-y2)^2)+@sqrt((x1-x3)^2+(y1-y3)^2)+@sqrt((x2-x3)^2+(y2-y3)^2);
(x1-5)^2+(y1-4)^2<=4;
(x2+5)^2+(y2+3)^2<=1;
(x3+1)^2+(y3-1)^2<=1;
@free(x1);
@free(y1);
@free(x2);
@free(y2);
@free(x3);
@free(y3);

结果是这个样子：

  Global optimal solution found.
  Objective value:                              18.41311
  Objective bound:                              18.41311
  Infeasibilities:                              0.000000
  Extended solver steps:                           43282
  Total solver iterations:                       4232645


                       Variable           Value        Reduced Cost
                             X1        3.361536            0.000000
                             X2       -4.180768          -0.2860243E-08
                             Y1        2.853075            0.000000
                             Y2       -2.426538            0.000000
                             X3      -0.2861699            0.000000
                             Y3       0.2996811            0.000000

                            Row    Slack or Surplus      Dual Price
                              1        18.41311           -1.000000
                              2        0.000000           0.5000000
                              3        0.000000            1.000000
                              4        0.000000            0.000000

可以看到，最小值是18.41311对应的坐标如上。

这个解起来挺复杂的，大概五六分钟的样子吧！嗯，像急性子的我可忍不了！让我们改进一下，用蒙特卡洛方法试一试。 解法如下： 1.在三个圆对应的矩形区域分别产生10万个模拟点（不高兴的话可以产生100万个或者更多）。 2.筛选出落在圆区域内的点。 3.分别从三个圆中选出一个点，计算它们组成的三角形的周长（总共10万个三角形）。 4.选出最小的周长对应的三角形，至此算法结束。 好了，Show me the code！ 对应的R语言程序如下（当然也可以用其他语言，只是R操作随机数更方便）：

rm(list = ls())
n<-100000#产生n个模拟点
x1<-runif(n,3,7)
y1<-runif(n,2,6)
x2<-runif(n,-6,-4)
y2<-runif(n,-4,-2)
x3<-runif(n,-2,0)
y3<-runif(n,0,2)
d<-data.frame(x1,y1,x2,y2,x3,y3)
r<-subset(d,(x1-5)^2+(y1-4)^2<=4&(x2+5)^2+(y2+3)^2<=1&(x3+1)^2+(y3-1)^2<=1)
#结果列表
rlist<-sqrt((r$x1-r$x2)^2+(r$y1-r$y2)^2)+sqrt((r$x1-r$x3)^2+(r$y1-r$y3)^2)+sqrt((r$x2-r$x3)^2+(r$y2-r$y3)^2)
rindex<-which.min(rlist)#最小值下标
cat("最小值=",rlist[rindex],"\n")
cat("x1=",r$x1[rindex],"\ny1=",r$y1[rindex],"\nx2=",r$x2[rindex],"\ny2=",r$y2[rindex],"\nx3=",r$x3[rindex],"\ny3=",r$y3[rindex])

这个的结果是这样的：

d=18.592
x1= 3.29677 
y1= 2.98348 
x2= -4.289648 
y2= -2.388063 
x3= -0.4919087 
y3= 0.2215853

差距不大吧，试着多运行几遍，每次都不一样，但最终收敛于最小值附近。 当然，上面的做法实际上大多数圆内的点都浪费了，因此，我们可以用参数方程的形式做得更好，具体的做法看代码：

rm(list = ls())
n<-100000#产生n个模拟点
theta<-runif(n,0,2*pi)#产生一个周期的角度
x1<-2*cos(theta)+5
y1<-2*sin(theta)+4
theta<-runif(n,0,2*pi)#产生一个周期的角度
x2<-cos(theta)-5
y2<-sin(theta)-3
theta<-runif(n,0,2*pi)#产生一个周期的角度
x3<-cos(theta)-1
y3<-sin(theta)+1
r<-data.frame(x1,x2,x3,y1,y2,y3)
rlist<-sqrt((r$x1-r$x2)^2+(r$y1-r$y2)^2)+sqrt((r$x1-r$x3)^2+(r$y1-r$y3)^2)+sqrt((r$x2-r$x3)^2+(r$y2-r$y3)^2)
rindex<-rindex<-which.min(rlist)#最小值下标
cat("最小值=",rlist[rindex],"\n")
cat("x1=",r$x1[rindex],"\ny1=",r$y1[rindex],"\nx2=",r$x2[rindex],"\ny2=",r$y2[rindex],"\nx3=",r$x3[rindex],"\ny3=",r$y3[rindex])

这样解出来的结果如下：

d=18.41324
x1= 3.366913 
y1= 2.845432 
x2= -4.183414 
y2= -2.422776 
x3= -0.4521829 
y3= 0.1634019

是不是跟实际的很接近了？这就是蒙特卡洛算法的奇特之处了。怎么样，感觉很easy对不对？