水資源管理分野の最適化問題：貯水池の放流スケジュール最適化

import numpy as np
from naginterfaces.library import opt

# パラメータ
I = np.array([3000, 4000, 15000, 50000, 30000, 12000, 8000, 5000, 4000, 3000, 3000, 4000, 3000, 5000])  # 各日の予測流入量（m³/日）
R_max = 15000           # 1日あたりの最大放流量（m³/日）
R_min = 1000            # 1日あたりの最小放流量（m³/日）
V_max = 2000000         # 貯水池の最大容量（m³）
V_min = 200000          # 貯水池の最小容量（m³）
V_max_opt = 1820000     # 貯水池の最適容量上限（m³）
V_min_opt = 1780000     # 貯水池の最適容量下限（m³）
P = 10                  # 1m³あたりの発電収益（円/m³）
C = 12                  # 貯水量が最適範囲を外れた場合のペナルティ（円/m³）

# 初期貯水量
V_start = (V_min_opt+V_max_opt)/2

# オプティマイザーの初期化
T = len(I)           # 計画期間（日数）
nvar = 3 * T
handle = opt.handle_init(nvar)

# 目的関数の係数（発電収益の最大化とペナルティの最小化）
cvec = [P] * T + [0] * T +  [-C] * T  # R, V, O の順番、Vは目的関数には含まれないのでゼロ
opt.handle_set_linobj(handle, cvec)

# 変数の上下限
bl = np.zeros(nvar)     # 下限の初期化
bu = np.zeros(nvar)     # 上限の初期化
bl[:T] = R_min          # R_t の下限
bu[:T] = R_max          # R_t の上限
bl[T:2*T] = V_min       # V_t の下限
bu[T:2*T] = V_max       # V_t の上限
bl[2*T:3*T] = 0         # O_t の下限
bu[2*T:3*T] = np.inf    # O_t の上限
opt.handle_set_simplebounds(handle, bl, bu)

# 初日の貯水量遷移制約
bl_1st_day = [V_start + I[0]]
bu_1st_day = [V_start + I[0]]
irowb_1st_day = [1, 1]
icolb_1st_day = [1, T+1]  # R1, V1
coeff_1st_day = [1, 1]    # R1 + V1
opt.handle_set_linconstr(handle, bl_1st_day, bu_1st_day, irowb_1st_day, icolb_1st_day, coeff_1st_day)

# 2日目以降の貯水量遷移制約
bl_2nd_day_on =I[1:]
bu_2nd_day_on =I[1:]
irowb_2nd_day_on = []
icolb_2nd_day_on = []
coeff_2nd_day_on = [1, -1, 1] * (T-1)
for t in range(1, T):
    bl_2nd_day_on[t-1] = I[t] 
    bu_2nd_day_on[t-1] = I[t] 
    irowb_2nd_day_on.extend([t,t,t])
    icolb_2nd_day_on.extend([t+1,T+t,T+t+1])  # Rt, Vt-1, Vt
opt.handle_set_linconstr(handle, bl_2nd_day_on, bu_2nd_day_on, irowb_2nd_day_on, icolb_2nd_day_on, coeff_2nd_day_on)

# 貯水適正量上限逸脱
bl_excess_upper = np.full(T, -V_max_opt)
bu_excess_upper = np.full(T, np.inf)
irowb_excess_upper = []
icolb_excess_upper = []
coeff_excess_upper = [1, -1] * T  # O_t - V_t
for t in range(T):
    irowb_excess_upper.extend([t+1, t+1])
    icolb_excess_upper.extend([2*T+t+1, T+t+1])  # O_t, V_t
opt.handle_set_linconstr(handle, bl_excess_upper, bu_excess_upper, irowb_excess_upper, icolb_excess_upper, coeff_excess_upper)

# 貯水適正量下限逸脱
bl_excess_lower = np.full(T, V_min_opt)
bu_excess_lower = np.full(T, np.inf)
irowb_excess_lower = []
icolb_excess_lower = []
coeff_excess_lower = [1, 1] * T
for t in range(T):
    irowb_excess_lower.extend([t+1, t+1])
    icolb_excess_lower.extend([2*T+t+1, T+t+1])  # O_t, V_t
opt.handle_set_linconstr(handle, bl_excess_lower, bu_excess_lower, irowb_excess_lower, icolb_excess_lower, coeff_excess_lower)

# ソルバーオプションの設定
opt.handle_opt_set(handle, 'Print Level = 0')
opt.handle_opt_set(handle, 'Task = Maximize')

# 初期値
x_init = [R_min]*T + [V_min]*T + [0]*T

# ソルバーの実行
optimal_solution, _, result_info, _ = opt.handle_solve_lp_ipm(handle,x_init)

#  結果を表示
R_opt = optimal_solution[:T]
V_opt = optimal_solution[T:2*T]
O_opt_upper = [max(0, v - V_max_opt) for v in V_opt]  # 上限超過量
O_opt_lower = [max(0, V_min_opt - v) for v in V_opt]  # 下限未満量

# 目的関数の値を取得
obj_value = result_info[1]

# 表示
print("\n最適解:")
print("日\t流入量\t\t放流量\t\t貯水量\t\t上限超過量\t下限未満量")
print(f"0\t-\t\t-\t\t{V_start:10.0f}\t-\t\t-")
for t in range(T):
    print(f"{t+1}\t{I[t]:<10}\t{R_opt[t]:10.0f}\t{V_opt[t]:10.0f}\t{O_opt_upper[t]:10.0f}\t{O_opt_lower[t]:10.0f}")

# 発電収益とペナルティー表示
revenue = sum(P * R_opt[t] for t in range(T))
penalty_upper = sum(C * O_opt_upper[t] for t in range(T))  # 上限超過ペナルティ
penalty_lower = sum(C * O_opt_lower[t] for t in range(T))  # 下限未満ペナルティ
penalty = penalty_upper + penalty_lower  # 総ペナルティ
print(f"-------------------------------")
print(f"総発電収益: {revenue:.0f}円")
print(f"上限超過ペナルティ: {penalty_upper:.0f}円")
print(f"下限未満ペナルティ: {penalty_lower:.0f}円")
print(f"総利益: {obj_value:.0f}円")