Home Energy Management¶

import numpy as np
    
num_days = 2
num_slots = num_days * 24
    
# Demands
# Temporal variation of predicted total demands (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
    
# Grid electricity
# Temporal variation of predicted electricity price in RTP (USD/kWh)
grid_electricity_price = np.array([
    [0.15, 0.14, 0.14, 0.13, 0.13, 0.14, 0.17, 0.21, 0.25, 0.28, 0.30, 0.32, # Day 1
     0.34, 0.36, 0.37, 0.39, 0.40, 0.42, 0.41, 0.38, 0.35, 0.32, 0.29, 0.25],
    [0.16, 0.15, 0.15, 0.14, 0.14, 0.15, 0.18, 0.22, 0.26, 0.29, 0.31, 0.33, # Day 2
     0.35, 0.37, 0.38, 0.40, 0.41, 0.43, 0.42, 0.39, 0.36, 0.33, 0.30, 0.27],
    ])  # fmt: skip
    
# Battery storage
battery_storage_capacity = 40  # Capacity of battery storage (kWh)
battery_storage_max_input_output = 20  # Maximum input/output of battery storage (kW)
    
# EV
battery_ev_capacity = 50  # Capacity of EV battery (kWh)
battery_ev_max_input = 50  # Maximum charge input of EV battery (kW)
# EV-chargeable hours (whether EV is parked at home (chargeable) at time t).
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 power
solar_panel_efficiency = 0.2  # Solar panel efficiency
solar_irradiance = 0.2  # Solar irradiance (kW/m^2)
panel_area = 20  # Panel area (m^2)
# Temporal variation of predicted sunshine hours (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):
    # Draw a vertical bar at every 24 hours
    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
    
    # Shade the time slots where EV is in use (not chargeable)
    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 in use",
                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 the problem setting
plot(
    {
        "demand (kW)": {"Total": np.ravel(demands)},
        "USD/kWh": {"grid electricity price": np.ravel(grid_electricity_price)},
        "supply (kW)": {"solar power": np.ravel(solar_power_max_supply)},
    },
)

import amplify
    
gen = amplify.VariableGenerator()
    
# Grid electricity
# q_ge corresponds to the grid electricity supply: 0, 1, .., 30 (kW)
q_ge = gen.array("Integer", shape=num_slots, bounds=(0, 30))
    
# Battery storage
# q_bs corresponds to the remaining battery charge: 0, 1, 2, ..., battery_storage_capacity (kWh)
q_bs = gen.array("Integer", shape=num_slots, bounds=(0, battery_storage_capacity))
# boundary conditions (50% of remaining charge at start and end of the optimization period)
q_bs[0] = int(0.5 * battery_storage_capacity)
q_bs[-1] = int(0.5 * battery_storage_capacity)
    
# Battery EV
# q_be corresponds to the remaining battery charge: 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)]
)
# Boundary conditions
q_be[0] = int(0.5 * battery_ev_capacity)
q_be[-1] = int(0.5 * battery_ev_capacity)
    
# EV battery is fully charged when leaving house
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
    
# Solar power
q_sp = gen.array("Binary", shape=num_slots)

supply_ge = q_ge  # kW
    
supply_bs = -(q_bs.roll(-1) - q_bs)  # kW, assuming time slot is 1 hour
    
supply_be = -(
    (q_be.roll(-1) - q_be) * np.roll(np.ravel(ev_parked_at_home), -1)
)  # kW, assuming time slot is 1 hour, EV battery output (when EV is in use) is not for household demand so forcing zero supply.
    
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 does not do positive supply, only charing
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()  # 1h
    - (supply_bs * supply_bs.roll(-2))[:-2].sum()  # 2h
    - (supply_bs * supply_bs.roll(-3))[:-3].sum()  # 3h
)

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:
    """A function that retries amplify.solve up to five times if no feasible solution is obtained."""
    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"  # If you use Amplify in a local environment, enter the Amplify API token.
    
# Constraint problem to obtain scaling factors
client.parameters.timeout = timedelta(milliseconds=5000)
# Obtain multiple solutions with the same objective value (there are multiple “optimal” solutions in Step 1, as no objective function is considered)
client.parameters.outputs.duplicate = True
model = amplify.Model(constraint_demand + constraint_bs + constraint_be)
result = retriable_solve(model, client)
    
# Calculate the objective function value for all feasible solutions, and consider
# the max value as the scaling factor for the corresponding objective function.
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}")

# Actual solution with full constraints and objectives
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(
    {
        "suuply/demand (kW)": {
            "total demand": np.ravel(demands),
            "total supply": (supply_ge + supply_bs + supply_be + supply_sp).evaluate(
                result.best.values
            ),
            "grid power": supply_ge.evaluate(result.best.values),
            "battery storage": supply_bs.evaluate(result.best.values),
            "battery EV": supply_be.evaluate(result.best.values),
            "solar power (kW)": supply_sp.evaluate(result.best.values),
        },
        "battery level (kWh)": {
            "battery storage": q_bs.evaluate(result.best.values),
            "battery EV": q_be.evaluate(result.best.values),
        },
    },
)