代码
#|echo: false
Sys.setenv(https_proxy = "http://127.0.0.1:7897", http_proxy = "http://127.0.0.1:7897")
2024年8月9日
线性优化(Linear Programming, LP)是一种数学优化方法,用于在特定约束条件下寻找线性目标函数的最优解。线性优化广泛应用于经济学、工程、交通、供应链管理、金融等众多领域。线性优化问题主要由2个部分组成。
目标函数
约束条件
考虑如下线性优化问题:
其中,目标函数是
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.
[1] "最优解"
[1] 0 6
[1] "目标函数值"
[1] 180
函数 OP()
定义一个优化问题,参数如下:
objective
:指定目标函数,用函数 L_objective()
表示线性优化中的目标函数,函数名中 L 表示 Linear(线性),包含数值型向量。constraints
:指定约束条件,用函数 L_constraint()
表示线性优化中的约束条件,函数名中 L 表示 Linear(线性),包含约束矩阵 >=
、<=
或 =
。types
:指定决策变量的类型,分三种情况, B
表示 0-1 变量,字母 B 是 binary 的意思,I
表示整型变量,字母 I 是 integer 的意思,C
表示数值型变量,字母 C 是 continuous 的意思。。maximum
:指定目标函数需要求极大还是极小,默认求极小,取值为逻辑值 TRUE
或 FALSE
。不同类型的目标函数和约束条件组合在一起可以构成非常丰富的优化问题。ROI 包支持的目标函数、约束条件及相应的代码见下表:
目标函数 | 代码 | 约束条件 | 代码 |
---|---|---|---|
线性函数 | L_objective() |
无约束 | 留空 |
二次函数 | Q_objective() |
箱式约束 | V_bound() |
非线性函数 | F_objective() |
线性约束 | L_constraint() |
二次约束 | Q_constraint() |
||
锥约束 | C_constraint() |
||
非线性约束 | F_constraint() |
旅行商问题(Travelling Salesman Problem,TSP)是一个经典的组合优化问题,它的目标是寻找一条最短路径,使得旅行者从一个起始城市出发,访问每一个城市一次且仅一次,最终返回到起点。这个问题在运筹学、计算机科学以及物流和运输等领域有着广泛的应用。
TSP包 是解决旅行商问题的最优工具包。该问题可被定义为:给定
可以将TSP问题等价于如下整数规划问题
从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 = "纬度")
通过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 。因采用启发式的随机优化算法,每次求解的结果也许会有所不同,建议采取不同的算法运行多次,选择最优的结果。
[1] 6133117
[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 = "纬度")
作为一个理性的投资者,希望回报最大而风险最小,给定投资和回报的约束条件下,选择风险最小的组合。一个简单的马科维茨投资组合优化问题如下:
其中,
下面基于股价数据介绍此组合优化问题。
首先利用 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
来求解组合优化问题。
[1] 0.4326 0.0872 0.2164 0.1190 0.0434 0.1014
[,1]
[1,] 0.7285183
求解出来的投资组合是43.26% 8.72% 21.64% 11.9% 4.34% 10.14%。
对于给定风险,最大化收益,可以建模如下
其中,目标函数中 Q_constraint()
来表示,这样线性约束和二次约束可以整合在一起,代码如下:
[,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()
quadprog
求解器不能求解该问题,尝试求解器 nloptr.slsqp
,不同股票同等看待,所以,权重的初始值都设置为
[1] 0.1615 0.3920 0.1763 0.1712 0.0991 0.0000
[,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.
[1] 0.1615 0.3920 0.1763 0.1712 0.0991 0.0000
[,1]
[1,] 0.1089005