import numpy as np
    
num_days = 2
num_slots = num_days * 24
    
# 需要
# 2日間の需要予測 (kW)
demands = np.array([
    #[0am,  1am,  2am,  3am,  4am,  5am,  6am,  7am,  8am,  9am, 10am, 11am,
    #12pm, 13pm, 14pm, 15pm, 16pm, 17pm, 18pm, 19pm, 20pm, 21pm, 22pm, 23pm]
    [1.28, 0.82, 1.96, 0.44, 0.88, 1.82, 5.70, 7.40, 2.64, 2.72, 1.52, 2.10, # Day 1
     1.20, 1.26, 4.00, 1.24, 2.36, 5.90, 5.92, 6.82, 8.86, 4.78, 1.08, 1.40],
    [1.56, 1.02, 0.40, 1.26, 1.72, 1.72, 4.52, 7.32, 5.66, 1.26, 3.68, 3.20, # Day 2
     3.54, 3.96, 1.28, 3.94, 2.98, 8.26, 5.5 , 6.46, 4.76, 6.32, 1.98, 0.58],
    ])  # fmt: skip
    
# 系統電力
# 2日間の RTP に基づく電力価格予測 (JPY/kWh)
grid_electricity_price = np.array([
    [24.0, 22.4, 22.4, 20.8, 20.8, 22.4, 27.2, 33.6, 40.0, 44.8, 48.0, 51.2, # Day 1
     54.4, 57.6, 59.2, 62.4, 64.0, 67.2, 65.6, 60.8, 56.0, 51.2, 46.4, 40.0],
    [25.6, 24. , 24. , 22.4, 22.4, 24. , 28.8, 35.2, 41.6, 46.4, 49.6, 52.8, # Day 2
     56.0, 59.2, 60.8, 64.0, 65.6, 68.8, 67.2, 62.4, 57.6, 52.8, 48.0, 43.2]
    ])  # fmt: skip
    
# 家庭用蓄電池
battery_storage_capacity = 40  # 最大容量 (kWh)
battery_storage_max_input_output = 20  # 最大入出力 (kW)
    
# EV バッテリー
battery_ev_capacity = 50  # 最大容量 (kWh)
battery_ev_max_input = 50  # 最大充電入力 (kW)
# EV 充電可能な時間帯の予測（時間スロット t に EV が家に駐車してある（充電可能）かどうか）.
ev_parked_at_home = np.array([
    [True, True, True, True, True, True, True, False, False, False, False, False,  # Day 1
     False, False, False, False, False, False, True, True, True, True, True, True, ],
    [True, True, True, True, True, True, True, False, False, False, False, False,  # Day 2
     False, False, False, False, False, False, True, True, True, True, True, True, ],
])  # fmt: skip
    
# 太陽光発電
solar_panel_efficiency = 0.2  # 太陽光パネル効率
solar_irradiance = 0.2  # 太陽放射照度 (kW/m^2)
panel_area = 20  # パネル面積 (m^2)
# 2日間の日照時間予測 (hours)
sunshine_hours = np.array([
    [0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.10, 0.30, 0.50, 0.70, 0.90, # Day 1
     1.00, 1.00, 0.90, 0.70, 0.50, 0.30, 0.10, 0.00, 0.00, 0.00, 0.00, 0.00],
    [0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.20, 0.40, 0.60, 0.80, 1.00, # Day 2
     0.90, 0.80, 0.70, 0.50, 0.30, 0.10, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00],
])  # hours  # fmt: skip
    
# 太陽光発電による最大供給可能量
solar_power_max_supply = (
    solar_panel_efficiency * solar_irradiance * panel_area * sunshine_hours
)  # kW
    
assert len(np.ravel(demands)) == num_slots
assert len(np.ravel(grid_electricity_price)) == num_slots
assert len(np.ravel(ev_parked_at_home)) == num_slots
assert len(np.ravel(sunshine_hours)) == num_slots

import plotly.graph_objects as go
from plotly.subplots import make_subplots
    
    
def plot(inp: dict[str, dict[str, np.ndarray]], show_ev_parked: bool = True):
    # 24時間ごとに目安の縦線を描画
    def daily_bar(fig: go.Figure) -> None:
        for i in range(num_days + 1):
            fig.add_vline(x=i * 24, line_dash="dot", line_color="#555555")
    
    def get_ev_usage_period(bools: list | np.ndarray) -> tuple[list, list]:
        start: list[float] = []
        end: list[float] = []
        if bools[0]:
            start.append(0)
        for i in range(len(bools) - 1):
            if not bools[i] and bools[i + 1]:
                start.append(i + 0.5)
            elif bools[i] and not bools[i + 1]:
                end.append(i + 0.5)
        if bools[-1]:
            end.append(len(bools) - 1)
        assert len(start) == len(end), f"length does not match {start=}, {end=}"
        return start, end
    
    # EV利用中（充電不可）の時間帯に影を描画
    def ev_usage(bools: list | np.ndarray, fig: go.Figure) -> None:
        start, end = get_ev_usage_period(bools)
        for s, e in zip(start, end):
            fig.add_vrect(
                x0=s,
                x1=e,
                fillcolor="black",
                opacity=0.1,
                annotation_text="EV利用中",
                annotation_position="top left",
            )
    
    num_plots = len(inp.keys())
    
    fig = make_subplots(rows=num_plots, cols=1)
    
    for i in range(num_plots):
        label, data = list(inp.items())[i]
        for k, v in data.items():
            fig.add_trace(
                go.Scatter(
                    x=list(range(len(v))),
                    y=v,
                    mode="lines",
                    name=f"{k}",
                ),
                row=1 + i,
                col=1,
            )
        fig.update_yaxes(title_text=label, row=1 + i, col=1)
    
    daily_bar(fig)
    if show_ev_parked:
        ev_usage([not elem for elem in np.ravel(ev_parked_at_home)], fig)
    fig.update_xaxes(title_text="hours", row=num_plots, col=1)
    fig.show()
    
    
# 問題設定をプロット
plot(
    {
        "需要 (kW)": {"総需要": np.ravel(demands)},
        "JPY/kWh": {"系統電力価格": np.ravel(grid_electricity_price)},
        "供給 (kW)": {"太陽光発電供給可能量": np.ravel(solar_power_max_supply)},
    },
)

import amplify
    
gen = amplify.VariableGenerator()
    
# 系統電力
# q_ge は系統電力供給: 0, 1, .., 30 (kW)
q_ge = gen.array("Integer", shape=num_slots, bounds=(0, 30))
    
# 家庭用蓄電池
# q_bs は充電残量: 0, 1, 2, ..., battery_storage_capacity (kWh)
q_bs = gen.array("Integer", shape=num_slots, bounds=(0, battery_storage_capacity))
# 境界条件 (最適化期間の開始・終了時に 50% の充電残量)
q_bs[0] = int(0.5 * battery_storage_capacity)
q_bs[-1] = int(0.5 * battery_storage_capacity)
    
# EV バッテリー
# q_be は充電残量: 0, 1, 2, ..., battery_ev_capacity (kWh)
q_be = gen.array("Integer", shape=num_slots, bounds=(0, battery_ev_capacity))
# pre-assign remaining batery charge to 10% if EV is not at home.
q_be = amplify.PolyArray(
    [q_be[i] if np.ravel(ev_parked_at_home)[i] else 0 for i in range(num_slots)]
)
# 境界条件 (最適化期間の開始・終了時に 50% の充電残量)
q_be[0] = int(0.5 * battery_ev_capacity)
q_be[-1] = int(0.5 * battery_ev_capacity)
    
# EV 利用開始時には満充電されている
for i in range(len(q_be) - 1):
    if np.ravel(ev_parked_at_home)[i] and not np.ravel(ev_parked_at_home)[i + 1]:
        q_be[i] = battery_ev_capacity
    
# 太陽光発電
q_sp = gen.array("Binary", shape=num_slots)

supply_ge = q_ge  # kW
    
supply_bs = -(q_bs.roll(-1) - q_bs)  # kW, 時間スロット幅が1時間であることを仮定
    
supply_be = -(
    (q_be.roll(-1) - q_be) * np.roll(np.ravel(ev_parked_at_home), -1)
)  # kW, 時間スロット幅が1時間であることを仮定
    
supply_sp = q_sp * np.ravel(solar_power_max_supply)  # kW

constraint_demand = amplify.ConstraintList(
    [
        amplify.greater_equal(
            (supply_ge + supply_bs + supply_be + supply_sp)[t],
            np.ravel(demands)[t],
            penalty_formulation="Relaxation",
        )
        for t in range(num_slots)
    ]
)
    
constraint_bs = amplify.ConstraintList(
    [
        amplify.less_equal(
            supply_bs[t] * supply_bs[t],
            battery_storage_max_input_output * battery_storage_max_input_output,
            penalty_formulation="Relaxation",
        )
        for t in range(num_slots)
    ]
)
    
# EV バッテリーは家庭電力需要に対して正の供給は行わない（充電のみ）。
constraint_be = amplify.ConstraintList(
    [amplify.clamp(supply_be[t], (-battery_ev_max_input, 0)) for t in range(num_slots)]
)

soft_constraint = (
    -(supply_bs * supply_bs.roll(-1))[:-1].sum()  # 1時間
    - (supply_bs * supply_bs.roll(-2))[:-2].sum()  # 2時間
    - (supply_bs * supply_bs.roll(-3))[:-3].sum()  # 3時間
)

objective_price = (supply_ge * np.ravel(grid_electricity_price)).sum()
    
supply_all = supply_ge + supply_bs + supply_be + supply_sp
objective_diff = (
    (np.ravel(demands) - supply_all) * (np.ravel(demands) - supply_all)
).sum()

from datetime import timedelta
    
    
def retriable_solve(
    model: amplify.Model, client: amplify.FixstarsClient
) -> amplify.Result:
    """ハード制約を全て満たす解が得られない場合は、5回まで求解をリトライする solve 関数"""
    for _ in range(5):
        result = amplify.solve(model, client)
        if len(result) > 0:
            break
        print("retrying amplify.solve...")
    if len(result) == 0:
        raise RuntimeError("no feasible solution found")
    return result
    
    
client = amplify.FixstarsClient()
# client.token = "xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx"  # ローカル環境等で使用する場合は、Fixstars Amplify AE のアクセストークンを入力してください。
    
# スケーリング係数を求めるための制約問題
client.parameters.timeout = timedelta(milliseconds=5000)
client.parameters.outputs.duplicate = True  # 同じ目的関数値の解を複数取得するオプション（目的関数を考慮しないステップ1では、“最適”解が複数個あるため）
model = amplify.Model(constraint_demand + constraint_bs + constraint_be)
result = retriable_solve(model, client)
    
# 得られた制約を満たす全ての解に対して目的関数値を計算し、その最大値を対応する目的関数のスケーリング係数とします。
num_sols = len(result.solutions)
scaling_diff = np.array(
    [objective_diff.evaluate(result.solutions[i].values) for i in range(num_sols)]
).max()
    
scaling_price = np.array(
    [objective_price.evaluate(result.solutions[i].values) for i in range(num_sols)]
).max()
    
scaling_soft = np.abs(
    np.array(
        [soft_constraint.evaluate(result.solutions[i].values) for i in range(num_sols)]
    )
).max()
    
print(f"{scaling_diff=:.10e}")
print(f"{scaling_price=:.10e}")
print(f"{scaling_soft=:.10e}")

# 全目的関数及び制約条件に基づく最適化問題
client.parameters.timeout = timedelta(milliseconds=10000)
model = amplify.Model(
    1000 * objective_diff / scaling_diff
    + 50 * objective_price / scaling_price
    + 0.1 * soft_constraint / scaling_soft,
    constraint_demand + constraint_bs + constraint_be,
)
    
result = retriable_solve(model, client)
    
print(f"objective_diff: {objective_diff.evaluate(result.best.values):.2e}")
print(f"objective_price: {objective_price.evaluate(result.best.values):.2e}")
print(f"soft_constraint: {soft_constraint.evaluate(result.best.values):.2e}")

plot(
    {
        "需要・供給 (kW)": {
            "総需要": np.ravel(demands),
            "総供給": (supply_ge + supply_bs + supply_be + supply_sp).evaluate(
                result.best.values
            ),
            "系統電力": supply_ge.evaluate(result.best.values),
            "家庭用蓄電池": supply_bs.evaluate(result.best.values),
            "EVバッテリー": supply_be.evaluate(result.best.values),
            "太陽光発電": supply_sp.evaluate(result.best.values),
        },
        "充電残量 (kWh)": {
            "家庭用蓄電池": q_bs.evaluate(result.best.values),
            "EVバッテリー": q_be.evaluate(result.best.values),
        },
    },
)

ホーム・エネルギー・マネジメントシステム¶

問題設定¶

定式化¶

決定変数¶

供給量¶

制約条件¶

ソフト制約¶

目的関数¶

求解の実行¶

ステップ 1. スケーリング係数を得るための求解¶

ステップ 2. 実際の求解¶

求解結果の可視化¶