强曰为道
与天地相似,故不违。知周乎万物,而道济天下,故不过。旁行而不流,乐天知命,故不忧.
文档目录

Julia 教程 / Julia 优化建模(JuMP.jl)

Julia 优化建模(JuMP.jl)

JuMP.jl 是 Julia 生态中的数学优化建模语言,支持线性规划、整数规划、非线性规划等多种优化问题。它提供简洁的代数建模语法,可以对接数十种优化求解器。


1. JuMP 概述

using JuMP, HiGHS

# 最简单的线性规划
# minimize 2x + 3y
# subject to: x + y >= 1, x >= 0, y >= 0
model = Model(HiGHS.Optimizer)

@variable(model, x >= 0)
@variable(model, y >= 0)
@objective(model, Min, 2x + 3y)
@constraint(model, x + y >= 1)

optimize!(model)

println("目标值: $(objective_value(model))")
println("x = $(value(x))")
println("y = $(value(y))")
println("求解状态: $(termination_status(model))")
核心宏说明
@variable(model, x)定义变量
@objective(model, sense, expr)定义目标函数
@constraint(model, expr)定义约束
@expression(model, expr)定义表达式(不创建约束)
@NLconstraint(model, expr)非线性约束
@NLobjective(model, sense, expr)非线性目标

2. 线性规划(LP)

using JuMP, HiGHS

# 生产计划问题
# 工厂生产两种产品,最大化利润
model = Model(HiGHS.Optimizer)

# 变量:产品产量
@variable(model, x1 >= 0)    # 产品 A
@variable(model, x2 >= 0)    # 产品 B

# 目标:最大化利润
@objective(model, Max, 40x1 + 30x2)

# 约束
@constraint(model, x1 + x2 <= 100)     # 原材料限制
@constraint(model, 2x1 + x2 <= 150)    # 工时限制
@constraint(model, x1 <= 60)           # 产品 A 需求上限
@constraint(model, x2 <= 80)           # 产品 B 需求上限

optimize!(model)

if termination_status(model) == OPTIMAL
    println("最优生产方案:")
    println("  产品 A: $(value(x1)) 件")
    println("  产品 B: $(value(x2)) 件")
    println("  最大利润: $(objective_value(model)) 元")

    # 灵敏度分析
    println("原材料影子价格: $(shadow_price(constraint_by_name(model, \"x1 + x2 <= 100\")))")
end

矩阵形式

using JuMP, HiGHS, LinearAlgebra

# 标准形式 LP: min c'x, s.t. Ax <= b, x >= 0
n, m = 5, 3
c = rand(n)
A = rand(m, n)
b = rand(m)

model = Model(HiGHS.Optimizer)
@variable(model, x[1:n] >= 0)
@objective(model, Min, dot(c, x))
@constraint(model, A * x .<= b)

optimize!(model)

3. 整数规划(IP)与混合整数规划(MIP)

using JuMP, HiGHS

# 背包问题:在容量限制下最大化价值
n = 10           # 物品数量
weights = rand(1:20, n)
values = rand(10:100, n)
capacity = 50

model = Model(HiGHS.Optimizer)

# 二进制变量:0-1 背包
@variable(model, x[1:n], Bin)

@objective(model, Max, dot(values, x))
@constraint(model, dot(weights, x) <= capacity)

optimize!(model)

selected = [i for i in 1:n if value(x[i]) > 0.5]
println("选中物品: $selected")
println("总价值: $(objective_value(model))")
println("总重量: $(sum(weights[selected]))")

MIP 示例:设施选址

using JuMP, HiGHS

# 设施选址问题
n_facilities = 5     # 候选设施数
n_customers = 10     # 客户数

# 成本参数
f = rand(50:200, n_facilities)    # 设施固定成本
c = rand(5:20, n_facilities, n_customers)  # 运输成本
d = rand(5:15, n_customers)       # 客户需求
cap = rand(50:100, n_facilities)  # 设施容量

model = Model(HiGHS.Optimizer)

# 决策变量
@variable(model, y[1:n_facilities], Bin)          # 是否开设设施
@variable(model, x[1:n_facilities, 1:n_customers] >= 0)  # 运输量

# 目标:最小化总成本
@objective(model, Min,
    sum(f[i] * y[i] for i in 1:n_facilities) +
    sum(c[i,j] * x[i,j] for i in 1:n_facilities, j in 1:n_customers)
)

# 容量约束
@constraint(model, [i in 1:n_facilities],
    sum(x[i,j] for j in 1:n_customers) <= cap[i] * y[i])

# 需求约束
@constraint(model, [j in 1:n_customers],
    sum(x[i,j] for i in 1:n_facilities) >= d[j])

optimize!(model)

opened = [i for i in 1:n_facilities if value(y[i]) > 0.5]
println("开设设施: $opened")
println("总成本: $(objective_value(model))")

4. 约束定义

using JuMP, HiGHS

model = Model(HiGHS.Optimizer)
@variable(model, x[1:5] >= 0)

# 单个约束
@constraint(model, sum(x) == 10)

# 范围约束
@constraint(model, 1 <= x[1] + x[2] <= 5)

# 索引约束族
@constraint(model, [i in 1:5], x[i] <= 3)

# 条件约束
@constraint(model, [i in 1:3, j in 4:5], x[i] + x[j] <= 5)

# 二次约束
@variable(model, y)
@constraint(model, x[1]^2 + x[2]^2 <= 1)    # 二阶锥约束

# 半正定约束
@variable(model, Z[1:3, 1:3], PSD)    # 半正定矩阵变量

# 命名约束
@constraint(model, budget, sum(x) <= 100)
shadow_price(budget)    # 获取影子价格

5. 目标函数

using JuMP, HiGHS

model = Model(HiGHS.Optimizer)
@variable(model, x[1:3] >= 0)

# 线性目标
@objective(model, Min, 2x[1] + 3x[2] + x[3])
@objective(model, Max, x[1] - x[2])

# 二次目标(二次规划 QP)
@objective(model, Min, x[1]^2 + 2x[2]^2 + x[3]^2 - x[1]*x[2])

# 无目标(可行性问题)
# 不调用 @objective 即可

# 修改目标
set_objective(model, Min, x[1] + x[2] + x[3])
set_objective_coefficient(model, x[1], 5.0)   # 修改 x[1] 的系数

6. 求解器选择

using JuMP

# 线性规划
model = Model(HiGHS.Optimizer)          # HiGHS(开源,推荐)
model = Model(GLPK.Optimizer)           # GLPK(开源)
model = Model(Clp.Optimizer)            # Clp(开源)
# model = Model(Gurobi.Optimizer)       # Gurobi(商业)
# model = Model(CPLEX.Optimizer)        # CPLEX(商业)

# 整数规划
model = Model(HiGHS.Optimizer)          # HiGHS 支持 MIP
model = Model(GLPK.Optimizer)           # GLPK 支持 MIP
# model = Model(Gurobi.Optimizer)       # Gurobi(最强 MIP 求解器)

# 二次规划(QP)
model = Model(Ipopt.Optimizer)          # Ipopt(非线性)
model = Model(OSQP.Optimizer)           # OSQP(凸 QP)

# 非线性规划
model = Model(Ipopt.Optimizer)          # Ipopt
model = Model(NLopt.Optimizer)          # NLopt

# 半正定规划
model = Model(COSMO.Optimizer)          # COSMO
求解器LPIP/MIPQPNLPSDP开源
HiGHS
GLPK
Ipopt
OSQP
COSMO
Gurobi
CPLEX

⚠️ 注意: 商业求解器(Gurobi、CPLEX)在 MIP 问题上性能远超开源求解器,但需要许可证。学术用户通常可以免费申请。

💡 提示: HiGHS 是目前最强的开源线性规划求解器,推荐作为默认选择。对于非线性问题,Ipopt 是首选。


7. 灵敏度分析

using JuMP, HiGHS

model = Model(HiGHS.Optimizer)
@variable(model, 0 <= x <= 10)
@variable(model, 0 <= y <= 10)
@objective(model, Max, 5x + 4y)
@constraint(model, con1, 6x + 4y <= 24)
@constraint(model, con2, x + 2y <= 6)

optimize!(model)

# 影子价格(对偶变量)
println("con1 影子价格: $(shadow_price(con1))")
println("con2 影子价格: $(shadow_price(con2))")

# 约束的右端项灵敏度范围
println("con1 RHS: $(normalized_rhs(con1))")
println("x 上界: $(upper_bound(x))")

# reduced cost(简约成本)
for (var, rc) in zip([x, y], [reduced_cost(x), reduced_cost(y)])
    println("$(name(var)) 的 reduced cost: $rc")
end

8. 多目标优化

using JuMP, HiGHS

# 加权和方法
model = Model(HiGHS.Optimizer)
@variable(model, x[1:3] >= 0)

# 两个冲突目标:最小化成本 vs 最大化质量
w = 0.6    # 权重
@objective(model, Min,
    w * (10x[1] + 15x[2] + 20x[3]) +            # 成本
    (1 - w) * -(5x[1] + 8x[2] + 12x[3])         # 质量(取负)
)

@constraint(model, sum(x) <= 100)
optimize!(model)

9. 建模技巧

9.1 线性化

# 非线性项 xy(其中 y 是二进制变量)的线性化
# xy 可以替换为 z,增加约束:
# z <= x, z <= M*y, z >= x - M*(1-y), z >= 0

using JuMP, HiGHS

M = 100    # big-M 常数

model = Model(HiGHS.Optimizer)
@variable(model, 0 <= x <= 50)
@variable(model, y, Bin)
@variable(model, z >= 0)

@constraint(model, z <= x)
@constraint(model, z <= M * y)
@constraint(model, z >= x - M * (1 - y))

@objective(model, Max, z)
optimize!(model)

9.2 分段线性函数

# 用 SOS2 约束建模分段线性函数
using JuMP, HiGHS

bp = [0.0, 10.0, 25.0, 50.0]
val = [0.0, 8.0, 20.0, 35.0]

model = Model(HiGHS.Optimizer)
@variable(model, x >= 0)
@variable(model, λ[1:length(bp)] >= 0)
@constraint(model, sum(λ) == 1)
@constraint(model, x == dot(λ, bp))
@constraint(model, λ in MOI.SOS2(collect(1:length(bp))))
@objective(model, Max, dot(λ, val))
optimize!(model)

10. 实际案例

10.1 运输问题

using JuMP, HiGHS

# 运输问题:最小化从仓库到客户的运输成本
n_supply = 3      # 仓库数
n_demand = 4      # 客户数

supply = [100, 150, 200]            # 各仓库供应量
demand = [80, 120, 100, 150]        # 各客户需求量
cost = [2 3 4 5;                    # 运输成本矩阵
        3 2 5 1;
        4 3 2 6]

# 平衡供需
@assert sum(supply) >= sum(demand)

model = Model(HiGHS.Optimizer)

@variable(model, x[1:n_supply, 1:n_demand] >= 0)

@objective(model, Min, sum(cost[i,j] * x[i,j] for i in 1:n_supply, j in 1:n_demand))

# 供应约束
@constraint(model, [i in 1:n_supply], sum(x[i,j] for j in 1:n_demand) <= supply[i])

# 需求约束
@constraint(model, [j in 1:n_demand], sum(x[i,j] for i in 1:n_supply) >= demand[j])

optimize!(model)

println("运输方案:")
for i in 1:n_supply
    for j in 1:n_demand
        if value(x[i,j]) > 0.01
            println("  仓库 $i → 客户 $j: $(value(x[i,j]))")
        end
    end
end
println("总运输成本: $(objective_value(model))")

10.2 排班问题

using JuMP, HiGHS

# 员工排班:7 天,每天需要不同人数
n_days = 7
n_employees = 15
required = [5, 4, 6, 7, 8, 6, 3]      # 每天需要最少人数
max_days = 5                            # 每人每周最多工作 5 天

model = Model(HiGHS.Optimizer)

@variable(model, x[1:n_employees, 1:n_days], Bin)

# 目标:最小化总工作天数
@objective(model, Min, sum(x))

# 每天人手约束
@constraint(model, [d in 1:n_days],
    sum(x[e, d] for e in 1:n_employees) >= required[d])

# 每人最多工作天数
@constraint(model, [e in 1:n_employees],
    sum(x[e, d] for d in 1:n_days) <= max_days)

# 至少连续休息 1 天
@constraint(model, [e in 1:n_employees, d in 1:n_days-1],
    x[e, d] + x[e, d+1] <= 1)    # 不能连续两天工作

optimize!(model)

println("排班表:")
for e in 1:n_employees
    schedule = [Int(value(x[e,d])) for d in 1:n_days]
    println("  员工 $e: $schedule")
end

10.3 投资组合优化

using JuMP, HiGHS, LinearAlgebra

# Markowitz 均值-方差投资组合
μ = [0.08, 0.12, 0.10, 0.06, 0.15]    # 预期收益率
Σ = [0.0064 0.0008 0.0011 0.0002 0.0018;   # 协方差矩阵
     0.0008 0.0144 0.0036 0.0012 0.0054;
     0.0011 0.0036 0.0100 0.0018 0.0045;
     0.0002 0.0012 0.0018 0.0049 0.0010;
     0.0018 0.0054 0.0045 0.0010 0.0225]

model = Model(HiGHS.Optimizer)
@variable(model, w[1:5] >= 0)              # 权重
@objective(model, Min, sum(Σ[i,j]*w[i]*w[j] for i in 1:5, j in 1:5))
@constraint(model, sum(w) == 1)
@constraint(model, dot(μ, w) >= 0.10)      # 最低收益 10%
optimize!(model)

println("最优权重: $(round.(value.(w)*100, digits=1))%")
println("预期收益: $(round(dot(μ, value.(w))*100, digits=2))%")

扩展阅读