R语言-运筹学非线性规划实例

Posted 2022-01-31 书槑

若目标函数或约束条件中含有自变量的非线性函数，则这样的规划问题就属于非线性规划。这里简单介绍4种情况下非线性规划的R语言实现。

目录

• 无约束规划
• 只含不等式约束的规划
• 系数可变的规划
• 约束条件含等式的规划

无约束规划

无约束规划中没有约束条件，只需要对目标函数取极值。R语言的nloptr包默认求解的是极小值，因此若要求最大值只需加个负号就行。

以这个目标函数为例，我们求它的极大值

library(nloptr) 
eval_f = function(x)
  return(list('objective'=-(12*x-3*x^4-2*x^6),   #输入目标函数
              'gradient'=-c(12-12*x^3-12*x^5)))  #一阶导要手动算一下

res0 = nloptr(x0 = c(0),  #给x赋一个初值，优化算法会不断迭代
              eval_f = eval_f,
              lb = c(-Inf),  #下界是负无穷
              ub = c(Inf),   #上界是正无穷
              opts = list('algorithm'='NLOPT_LD_LBFGS'))
res0$solution

输出的结果是x的值，为局部最优

但由于这是一个凸规划，任何局部最优解也是其全局最优解。

只含不等式约束的规划

在用R语言求解时，需要加上约束条件的方程，默认的不等式约束条件是小于等于0。还要给出雅各比矩阵，即变量的一阶偏导矩阵。

eval_f = function(x)
  return(list('objective'=-(3*x[1]+5*x[2]),
              'gradient'=c(-3,-5)))


eval_g = function(x)
  return(list('constraints' = c(8*x[1]-x[1]^2+14*x[2]-x[2]^2-49,x[1]-4),
              'jacobian' = rbind(c(8-2*x[1],14-2*x[2]),c(1,0))))


res1 = nloptr(x0 = c(0,0),
              eval_f = eval_f,
              lb = c(0,0),
              ub = c(4,Inf),
              eval_g_ineq = eval_g,
              opts = list('algorithm'='NLOPT_LD_MMA'))
res1$solution  

这样求出来的结果也是局部最优，若改变初始值，结果可能就不一样了。

系数可变的规划

在这种情况下，约束条件中的系数是可变的

题目如下

# min sqrt( x2 )
# s.t. x2 >= 0
#      x2 >= ( a1*x1 + b1 )^3
#      x2 >= ( a2*x1 + b2 )^3
# where
# a1 = 2, b1 = 0, a2 = -1, b2 = 1

求解过程如下

a = c(2,-1)
b = c(0,1)
eval_f = function(x,a,b)
  return(list('objective'=sqrt(x[2]),
              'gradient'=c(0,1/(2*sqrt(x[2])))))

eval_g = function(x,a,b)
  return(list('constraints' = c((a[1]*x[1]+b[1])^3-x[2],
                               (a[2]*x[1]+b[2])^3-x[2]),
              'jacobian' = rbind(c(3*a[1]*(a[1]*x[1]+b[1])^2,-1),
                                 c(3*a[2]*(a[2]*x[1]+b[2])^2,-1))))

res2 = nloptr(x0 = c(1.234,5.678),
              eval_f = eval_f,
              lb = c(-Inf,0),
              ub = c(Inf,Inf),
              eval_g_ineq = eval_g,
              opts = list('algorithm'='NLOPT_LD_MMA'),
              a = a, b = b)
res2$solution

约束条件含等式的规划

等式约束需单独用一个函数写进去，同样要写出相应的雅各比矩阵

题目如下

#min x1^2 + x2^2 + 8
#s.t. x1^2 - x2 >= 0
#   -x1 - x2^2 + 2 = 0
#   x1, x2 >= 0

求解过程如下

eval_f = function(x)
  return(list('objective'=x[1]^2 + x[2]^2 + 8,
              'gradient'=c(2*x[1],2*x[2])))

eval_g = function(x)
  return(list('constraints' = c(x[2]-x[1]^2),
              'jacobian' = rbind(c(-2*x[1],1))))

eval_h = function(x)
  return(list('constraints' = c(-x[1]-x[2]^2+2),
              'jacobian' = rbind(c(-1,-2*x[2]))))

res3 = nloptr(x0 = c(2,0),
              eval_f = eval_f,
              lb = c(0,0),
              ub = c(Inf,Inf),
              eval_g_ineq = eval_g,
              eval_g_eq = eval_h,
              opts = list('algorithm'='NLOPT_LD_SLSQP'))
res3$solution