优化问题

一些运筹优化相关问题的汇总(旅行商问题、投资组合问题)

#|echo: false
Sys.setenv(https_proxy = "http://127.0.0.1:7897", http_proxy = "http://127.0.0.1:7897")

线性优化

线性优化（Linear Programming, LP）是一种数学优化方法，用于在特定约束条件下寻找线性目标函数的最优解。线性优化广泛应用于经济学、工程、交通、供应链管理、金融等众多领域。线性优化问题主要由2个部分组成。

• 目标函数 $$min/maxz=c_1{x}_1+c_2{x}_2+...+c_n{x}_n$$

• 约束条件 $$\begin{array}{rl}\text{s.t.}& \{\begin{array}{c}{a}_{11}{x}_1+a_{12}{x}_2+a_{1n}{x}_n\le b_1\\ a_{21}{x}_1+a_{22}{x}_2+a_{2n}{x}_n\le b_2\\ ...\\ a_{n1}{x}_1+a_{n2}{x}_2+a_{nn}{x}_n\le b_n\end{array}\end{array}$$

考虑如下线性优化问题：

$$\begin{array}{rl}\underset{\mathit{x}}{max}& 40x_1+30x_2\\ \text{s.t.}& \{\begin{array}{c}2x_1+1x_2\le 8\\ 3x_1+2x_2\le 12\\ 0\le x_1,0\le x_2\end{array}\end{array}$$

其中，目标函数是 $40x_1+30x_2$，$max$ 表示求目标函数的最小值$\text{s.t.}$ 是 subject to 的缩写，专指约束条件。 用 ROI 包提供的一套使用语法表示该线性优化问题，代码如下：

library(ROI)
# 定义优化问题
op <- OP(
  objective = L_objective(L = c(40, 30)),
  constraints = L_constraint(
    L = matrix(c(
      2, 1,
      3, 2,
      1, 0,
      0, 1
    ), ncol = 2, byrow = TRUE),
    dir = c("<=", "<=", ">=", ">="),
    rhs = c(8, 12, 0, 0)
  ),
  types = c("C", "C"),
  maximum = T
)
# 优化问题描述
op

ROI Optimization Problem:

Maximize a linear objective function of length 2 with
- 2 continuous objective variables,

subject to
- 4 constraints of type linear.
- 0 lower and 0 upper non-standard variable bounds.

# 求解优化问题
res <- ROI_solve(op, solver = "glpk")
print("最优解")

[1] "最优解"

res$solution

[1] 0 6

print("目标函数值")

[1] "目标函数值"

res$objval

[1] 180

函数 OP() 定义一个优化问题，参数如下：

• objective ：指定目标函数，用函数 L_objective() 表示线性优化中的目标函数，函数名中 L 表示 Linear（线性），包含数值型向量。
• constraints ：指定约束条件，用函数 L_constraint() 表示线性优化中的约束条件，函数名中 L 表示 Linear（线性），包含约束矩阵 $A$ ，约束分量的方向可为 >= 、<= 或 = 。
• types ：指定决策变量的类型，分三种情况， B 表示 0-1 变量，字母 B 是 binary 的意思，I 表示整型变量，字母 I 是 integer 的意思，C 表示数值型变量，字母 C 是 continuous 的意思。。
• maximum ：指定目标函数需要求极大还是极小，默认求极小，取值为逻辑值 TRUE 或 FALSE。

不同类型的目标函数和约束条件组合在一起可以构成非常丰富的优化问题。ROI 包支持的目标函数、约束条件及相应的代码见下表:

ROI 包可以表示的目标函数和约束条件

目标函数	代码	约束条件	代码
线性函数	L_objective()	无约束	留空
二次函数	Q_objective()	箱式约束	V_bound()
非线性函数	F_objective()	线性约束	L_constraint()
		二次约束	Q_constraint()
		锥约束	C_constraint()
		非线性约束	F_constraint()

旅行商问题

旅行商问题（Travelling Salesman Problem，TSP）是一个经典的组合优化问题，它的目标是寻找一条最短路径，使得旅行者从一个起始城市出发，访问每一个城市一次且仅一次，最终返回到起点。这个问题在运筹学、计算机科学以及物流和运输等领域有着广泛的应用。

TSP包 是解决旅行商问题的最优工具包。该问题可被定义为：给定 $n$个城市间的距离，用矩阵$D$表示，其中元素$d_{ij}$表示城市$i$到城市$j$ 之间的距离，且对角元素$d_{ii}=0$（即城市到自身的距离为零），其中 $i,j=1,2,\cdots ,n$ 。一条旅行路线可以表示为循环排 $\{1,2,\dots ,n\}$的某种排列$\pi$，其中$\pi (i)$表示在城市$i$之后的城市。旅行商问题的目标是找到一个排列$\pi$，使得对应的旅行线路总距离最短

$$\sum _{i=1}^n{d}_{i\pi (i)}$$

可以将TSP问题等价于如下整数规划问题

$$\begin{array}{rl}min& \sum _{i=1}^n\sum _{j=1}^n{d}_{ij}{x}_{ij}\\ \text{s.t.}& \sum _{i=1}^n{x}_{ij}=1,j=1,2,\dots ,n,\\ & \sum _{j=1}^n{x}_{ij}=1,i=1,2,\dots ,n,\\ & x_{ij}=0\text{or}1\end{array}$$

示例

从maps包可以获取中国人口最多的10个城市的坐标，从北京出发，最后回到北京，如何规划线路才能使得总行程最短？行程最短的路径是什么？10 个城市的分布如 图 1 所示。

# 获取坐标数据
library(dplyr)
city <- maps::world.cities %>%
  filter(country.etc == "China") %>%
  slice_max(n = 10, order_by = pop)


library(sf)
# crs = 4326，表示指定采用 WGS 84 地理坐标系
city_latlong <- st_as_sf(city,
  coords = c("long", "lat"), crs = 4326
)

library(ggplot2)
ggplot() +
  geom_sf_label(
    data = city_latlong, aes(label = name),
    fun.geometry = sf::st_centroid
  ) +
  geom_sf(data = city_latlong, color = "red") +
  theme_minimal() +
  labs(x = "经度", y = "纬度")

图 1: 10 个城市的分布图

通过sf包的st_distance可以计算城市坐标之间的直线距离

# 计算距离矩阵
distance_matrix <- st_distance(city_latlong)
distance_matrix <- matrix(distance_matrix, nrow = length(city_latlong$name))
rownames(distance_matrix) <- city_latlong$name
colnames(distance_matrix) <- city_latlong$name

library(TSP)
D_tsp <- as.TSP(distance_matrix)

tour_sol <- solve_TSP(x = D_tsp, method = "nearest_insertion", start = 2)
tour_sol

object of class 'TOUR' 
result of method 'nearest_insertion' for 10 cities
tour length: 6133117 

途经 10 个城市的最短路程为 6133117 m 。因采用启发式的随机优化算法，每次求解的结果也许会有所不同，建议采取不同的算法运行多次，选择最优的结果。

# 旅行最短路程
tour_length(tour_sol)

[1] 6133117

# 旅行线路方案
labels(D_tsp)[as.integer(tour_sol)]

 [1] "Beijing"   "Harbin"    "Shenyang"  "Shanghai"  "Nanjing"   "Wuhan"    
 [7] "Chongqing" "Chengdu"   "Xian"      "Tianjin"  

求解结果对应的旅行方案，如 图 2 所示，依次走过的城市是：北京、天津、西安、成都、重庆、武汉、南京、上海、沈阳、哈尔滨。

city_tour <- st_cast(st_combine(st_geometry(city_latlong[as.integer(tour_sol), ])), "POLYGON")
ggplot() +
  geom_sf_label(
    data = city_latlong, aes(label = name),
    fun.geometry = sf::st_centroid
  ) +
  geom_sf(data = city_latlong, color = "steelblue") +
  geom_sf(data = city_latlong %>% filter(name == "Beijing"), color = "red", size = 3) +
  geom_sf(data = city_tour, fill = NA, color = "black") +
  theme_minimal() +
  labs(x = "经度", y = "纬度")

图 2: 10 个城市的路线图

投资组合问题

作为一个理性的投资者，希望回报最大而风险最小，给定投资和回报的约束条件下，选择风险最小的组合。一个简单的马科维茨投资组合优化问题如下：

$$\begin{array}{rl}\underset{\mathit{w}}{min}& {\mathit{w}}^{\mathrm{\top }}\hat{\mathrm{\Sigma }}\mathit{w}\\ \text{s.t.}& A{\mathit{w}}^{\mathrm{\top }}\le \mathit{b}\end{array}$$

其中，$\mathit{w}$ 是权重向量，每个分量代表对投资对象的投资比例，$\hat{\mathrm{\Sigma }}$ 是关于投资对象的协方差矩阵，约束条件中包含两个部分，一个是权重之和为 1，一个是投资组合的收益率达到预期值。

示例

下面基于股价数据介绍此组合优化问题。

首先利用 tidyquant 包获取2024年后的历史股价数据，再根据 期间的股票调整价，使用tq_transmute计算各支股票天粒度的收益率。收益率可以看作一个随机变量，收益率的波动变化，即随机变量的方差，可以看作风险。

对于给定收益率，最小化风险情况，可以建模如下

# 加载需要的库
library(ROI) # 优化求解接口
library(ROI.plugin.glpk) # GLPK（GNU线性规划工具包）插件，用于线性规划
library(ROI.plugin.nloptr) # NLOPTR 插件，用于非线性优化
library(ROI.plugin.scs) # SCS 插件，用于稀疏约束最优化
library(ROI.plugin.quadprog) # Quadprog 插件，用于二次规划
library(lattice) # 用于数据可视化
library(tidyquant) # 用于金融数据的获取和处理
library(tibble) # 用于数据框操作

# 从金融数据源获取指定股票的价格数据
tb_stocks <- tq_get(c("AMZN", "AAPL", "INTC", "IBM", "000001.ss", "600036.ss"),
  get = "stock.prices", from = "2024-01-01", to = "2024-06-01"
)

# 计算每日收益率
tb_stocks_return <- tb_stocks %>%
  group_by(symbol) %>% # 按股票符号分组
  tq_transmute(adjusted, periodReturn, period = "daily", col_rename = "return") %>%
  pivot_table(.columns = ~symbol, .values = ~return, .rows = ~date) %>% # 透视表，将收益率按日期和股票符号排列
  column_to_rownames("date") %>% # 将日期列转换为行名
  na.omit() # 删除含有 NA 的行

DD <- 100 * as.matrix(tb_stocks_return)

# 计算平均收益率
r <- mean(DD, na.rm = TRUE) # 计算每日收益率的平均值
r # 输出平均收益率

[1] -0.004641027

# 定义目标函数，使用协方差矩阵作为目标
foo <- Q_objective(Q = cov(DD), L = rep(0, ncol(DD))) # 二次规划目标函数

# 设定投资组合的总投资约束
full_invest <- L_constraint(rep(1, ncol(DD)), "==", 1) # 投资比例总和等于1

# 设定回报约束，确保回报达到平均水平
target_return <- L_constraint(apply(DD, 2, mean), "==", r) # 回报水平约束

# 创建优化问题
op <- OP(objective = foo, constraints = rbind(full_invest, target_return)) # 组合目标函数和约束条件


op

ROI Optimization Problem:

Minimize a quadratic objective function of length 6 with
- 6 continuous objective variables,

subject to
- 2 constraints of type linear.
- 0 lower and 0 upper non-standard variable bounds.

求解器 nloptr.slsqp 需要给初值和等式约束的梯度，而求解器 quadprog 不需要给初值。下面使用 quadprog 来求解组合优化问题。

library(ROI.plugin.quadprog)
sol <- ROI_solve(op, solver = "quadprog")
# 最优解：投资组合
w <- sol$solution
# 保留 4 位小数
round(w, 4)

[1] 0.4326 0.0872 0.2164 0.1190 0.0434 0.1014

# 目标函数值：投资风险
sqrt(t(w) %*% cov(DD) %*% w)

          [,1]
[1,] 0.7285183

求解出来的投资组合是43.26% 8.72% 21.64% 11.9% 4.34% 10.14%。

对于给定风险，最大化收益，可以建模如下

$$\begin{array}{rl}\underset{\mathit{w}}{max}& {\mathit{w}}^{\mathrm{\top }}\hat{\mathit{\mu }}\\ \text{s.t.}& A\mathit{w}\le \mathit{b}\\ & {\mathit{w}}^{\mathrm{\top }}\hat{\mathrm{\Sigma }}\mathit{w}\le \sigma \end{array}$$

其中，目标函数中 $\hat{\mathit{\mu }}$ 表示根据历史数据获得的投资对象的收益率，约束条件中 $\sigma$ 表示投资者可以接受的投资风险，其他符号的含义同前。在给定风险约束 $\sigma$ 下，求取回报最大的组合。线性约束也可以用函数 Q_constraint() 来表示，这样线性约束和二次约束可以整合在一起，代码如下：

# 风险阈值
sigma <- sqrt(t(w) %*% cov(DD) %*% w)
sigma

          [,1]
[1,] 0.7285183

# 全 0 矩阵
zero_mat <- diag(x = rep(0, ncol(DD)))
# 目标函数
foo <- Q_objective(Q = zero_mat, L = colMeans(DD))
# 线性和二次约束
maxret_constr <- Q_constraint(
  Q = list(cov(DD), NULL),
  L = rbind(
    rep(0, ncol(DD)),
    rep(1, ncol(DD))
  ),
  dir = c("<=", "=="), rhs = c(1 / 2 * sigma^2, 1)
)
# 目标规划
op <- OP(objective = foo, constraints = maxret_constr, maximum = TRUE)
op

ROI Optimization Problem:

Maximize a quadratic objective function of length 6 with
- 6 continuous objective variables,

subject to
- 2 constraints of type quadratic.
- 0 lower and 0 upper non-standard variable bounds.

ROI_applicable_solvers()

函数 ROI_applicable_solvers() 识别规划问题类型，给出可求解此规划问题的求解器。

ROI_applicable_solvers(op)

[1] "nloptr.cobyla" "nloptr.mma"    "nloptr.auglag" "nloptr.isres" 
[5] "nloptr.slsqp" 

quadprog 求解器不能求解该问题，尝试求解器 nloptr.slsqp ，不同股票同等看待，所以，权重的初始值都设置为 $\frac{1}{6}$ 。

# 求解规划问题
nlp <- ROI_solve(op, solver = "nloptr.slsqp", start = rep(1 / 6, 6))
# 投资组合
w <- nlp$solution
# 保留 4 位小数
round(w, 4)

[1] 0.1615 0.3920 0.1763 0.1712 0.0991 0.0000

# 投资组合的预期收益
w %*% colMeans(DD)

          [,1]
[1,] 0.1089005

包含一些非线性的等式或不等式约束,可以用函数 F_constraint() 表示的代码如下

# x 是一个表示权重的列向量
# 等式约束
# 权重之和为 1 的约束
heq <- function(x) {
  sum(x)
}
# 等式约束的雅可比
heq.jac <- function(x) {
  rep(1, length(x))
}
# 不等式约束
# 二次的风险约束
hin <- function(x) {
  1 / 2 * t(x) %*% cov(DD) %*% x
}
# 不等式约束的雅可比
hin.jac <- function(x) {
  cov(DD) %*% x
}
# 目标规划
op <- OP(
  objective = L_objective(L = colMeans(DD)), # 6 个目标变量
  constraints = F_constraint(
    # 等式和不等式约束
    F = list(heq = heq, hin = hin),
    dir = c("==", "<="),
    rhs = c(1, 1 / 2 * sigma^2),
    # 等式和不等式约束的雅可比
    J = list(heq.jac = heq.jac, hin.jac = hin.jac)
  ),
  # 目标变量的取值范围
  bounds = V_bound(ld = 0, ud = 1, nobj = 6L),
  maximum = TRUE # 最大回报
)
op

ROI Optimization Problem:

Maximize a linear objective function of length 6 with
- 6 continuous objective variables,

subject to
- 2 constraints of type nonlinear.
- 0 lower and 6 upper non-standard variable bounds.

# 求解规划问题
nlp <- ROI_solve(op, solver = "nloptr.slsqp", start = rep(1 / 6, 6))
# 投资组合
w <- nlp$solution
round(w, 4)

[1] 0.1615 0.3920 0.1763 0.1712 0.0991 0.0000

# 投资组合的预期收益
w %*% colMeans(DD)

          [,1]
[1,] 0.1089005