Python Smart Grid Simulation — Deep Dive

Build production-grade smart grid simulations in Python using power flow solvers, optimization models, and reinforcement learning for grid control.

Technical foundation

Smart grid simulation encompasses power flow analysis, optimal dispatch, network planning, and automated control. Python serves as the integration layer connecting physics-based solvers with optimization algorithms and machine learning controllers. This deep dive covers the computational methods, code patterns, and architectural decisions for building realistic grid simulations.

Network modeling with pandapower

pandapower provides a pandas-based interface for creating and solving power system models:

import pandapower as pp
import pandapower.networks as nw

# Create a simple distribution network
net = pp.create_empty_network()

# Add buses (connection points)
bus_hv = pp.create_bus(net, vn_kv=110, name="HV Bus")
bus_mv = pp.create_bus(net, vn_kv=20, name="MV Bus")
bus_lv1 = pp.create_bus(net, vn_kv=0.4, name="LV Bus 1")
bus_lv2 = pp.create_bus(net, vn_kv=0.4, name="LV Bus 2")

# Add transformer (HV to MV)
pp.create_transformer(net, bus_hv, bus_mv, std_type="63 MVA 110/20 kV")

# Add MV to LV transformers
pp.create_transformer(net, bus_mv, bus_lv1, std_type="0.4 MVA 20/0.4 kV")
pp.create_transformer(net, bus_mv, bus_lv2, std_type="0.4 MVA 20/0.4 kV")

# Add lines
pp.create_line(net, bus_mv, bus_mv, length_km=2.0, std_type="NA2XS2Y 1x185 RM/25 12/20 kV")

# Add external grid connection (slack bus)
pp.create_ext_grid(net, bus_hv, vm_pu=1.02)

# Add loads
pp.create_load(net, bus_lv1, p_mw=0.2, q_mvar=0.05, name="Residential Area 1")
pp.create_load(net, bus_lv2, p_mw=0.15, q_mvar=0.03, name="Residential Area 2")

# Add distributed solar generation
pp.create_sgen(net, bus_lv1, p_mw=0.08, q_mvar=0, name="Rooftop Solar 1")
pp.create_sgen(net, bus_lv2, p_mw=0.12, q_mvar=0, name="Rooftop Solar 2")

# Run power flow
pp.runpp(net, algorithm="nr")  # Newton-Raphson

# Check results
print(net.res_bus)   # Voltage magnitude and angle per bus
print(net.res_line)  # Line loading percentage
print(net.res_trafo) # Transformer loading

Key checks after power flow:

  • Voltage limits: All buses within 0.95–1.05 p.u. (per unit of nominal).
  • Line loading: Below 100% thermal rating (80% for N-1 contingency).
  • Transformer loading: Below rated capacity with margin for contingency.

Time-series simulation

Annual simulations require driving the network model with time-varying load and generation profiles:

import numpy as np
import pandas as pd
import pandapower.timeseries as ts
from pandapower.timeseries.data_sources.frame_data import DFData

# Create hourly profiles for a year (8760 hours)
hours = 8760
time_index = pd.date_range("2025-01-01", periods=hours, freq="1h")

# Load profile: base load + temperature-dependent component
np.random.seed(42)
base_load = 0.15  # MW
temp_profile = 15 + 10 * np.sin(2 * np.pi * np.arange(hours) / 8760 - np.pi/2)
load_profile = base_load + 0.05 * np.maximum(temp_profile - 22, 0)  # Cooling
load_profile += np.random.normal(0, 0.01, hours)  # Noise

# Solar profile: clear-sky with cloud variability
solar_max = 0.12  # MW peak
hour_of_day = np.arange(hours) % 24
solar_clear = solar_max * np.maximum(np.sin(np.pi * (hour_of_day - 6) / 12), 0)
cloud_factor = np.random.uniform(0.3, 1.0, hours)
solar_profile = solar_clear * cloud_factor

# Create data sources
load_data = DFData(pd.DataFrame({"p_mw": load_profile}, index=time_index))
solar_data = DFData(pd.DataFrame({"p_mw": solar_profile}, index=time_index))

# Create controllers for time-varying elements
ts.create_output_writer(net, time_steps=hours)
# Controller binds data source to network element
const_load = ts.ConstControl(net, element="load", variable="p_mw",
                              element_index=[0], data_source=load_data,
                              profile_name="p_mw")
const_solar = ts.ConstControl(net, element="sgen", variable="p_mw",
                               element_index=[0], data_source=solar_data,
                               profile_name="p_mw")

# Run time series
ts.run_timeseries(net, time_steps=range(hours))

Optimal power flow and economic dispatch

Beyond checking physical feasibility, optimal power flow (OPF) minimizes generation cost while respecting network constraints:

import pyomo.environ as pyo

def economic_dispatch(generators: list[dict], demand: float) -> dict:
    """Simple economic dispatch using quadratic cost curves."""
    model = pyo.ConcreteModel()

    model.G = pyo.Set(initialize=range(len(generators)))
    model.p = pyo.Var(model.G, within=pyo.NonNegativeReals)

    # Objective: minimize total generation cost
    def cost_rule(m):
        return sum(
            generators[g]["a"] * m.p[g]**2
            + generators[g]["b"] * m.p[g]
            + generators[g]["c"]
            for g in m.G
        )
    model.cost = pyo.Objective(rule=cost_rule, sense=pyo.minimize)

    # Power balance constraint
    model.balance = pyo.Constraint(
        expr=sum(model.p[g] for g in model.G) == demand
    )

    # Generator limits
    def pmin_rule(m, g):
        return m.p[g] >= generators[g]["p_min"]
    def pmax_rule(m, g):
        return m.p[g] <= generators[g]["p_max"]
    model.pmin = pyo.Constraint(model.G, rule=pmin_rule)
    model.pmax = pyo.Constraint(model.G, rule=pmax_rule)

    solver = pyo.SolverFactory("glpk")
    solver.solve(model)

    return {g: pyo.value(model.p[g]) for g in model.G}

For smart grids, the dispatch problem extends to include battery scheduling, demand response activation, and renewable curtailment.

Battery storage optimization

Battery dispatch involves deciding when to charge and discharge over a planning horizon:

def optimize_battery_schedule(
    prices: np.ndarray,      # Electricity price per hour
    solar: np.ndarray,        # Solar generation per hour
    load: np.ndarray,         # Load demand per hour
    battery_capacity: float,  # kWh
    max_power: float,         # kW charge/discharge rate
    efficiency: float = 0.92, # Round-trip efficiency
    hours: int = 24,
) -> dict:
    """Optimize battery charge/discharge to minimize electricity cost."""
    model = pyo.ConcreteModel()
    model.T = pyo.RangeSet(0, hours - 1)

    model.charge = pyo.Var(model.T, bounds=(0, max_power))
    model.discharge = pyo.Var(model.T, bounds=(0, max_power))
    model.soc = pyo.Var(model.T, bounds=(0.1 * battery_capacity, 0.9 * battery_capacity))

    # Grid import = load - solar - discharge + charge
    def grid_import(m, t):
        return load[t] - solar[t] - m.discharge[t] + m.charge[t]

    # Minimize cost of grid imports
    model.obj = pyo.Objective(
        expr=sum(prices[t] * grid_import(model, t) for t in model.T),
        sense=pyo.minimize,
    )

    # SOC dynamics
    def soc_rule(m, t):
        if t == 0:
            return m.soc[t] == battery_capacity * 0.5  # Start at 50%
        return m.soc[t] == (
            m.soc[t-1]
            + m.charge[t] * efficiency**0.5
            - m.discharge[t] / efficiency**0.5
        )
    model.soc_dynamics = pyo.Constraint(model.T, rule=soc_rule)

    solver = pyo.SolverFactory("glpk")
    solver.solve(model)

    return {
        "charge": [pyo.value(model.charge[t]) for t in model.T],
        "discharge": [pyo.value(model.discharge[t]) for t in model.T],
        "soc": [pyo.value(model.soc[t]) for t in model.T],
        "cost": pyo.value(model.obj),
    }

Reinforcement learning for grid control

Grid2Op provides a gym-compatible environment for training RL agents to operate power grids:

import grid2op
from grid2op.Agent import DoNothingAgent

# Create environment
env = grid2op.make("l2rpn_case14_sandbox")

# Simple agent that reconnects disconnected lines
class ReconnectAgent(grid2op.Agent.BaseAgent):
    def act(self, observation, reward, done=False):
        action = self.action_space({})
        # Find disconnected lines and try to reconnect
        disconnected = ~observation.line_status
        if disconnected.any():
            line_id = disconnected.nonzero()[0][0]
            action = self.action_space({"set_line_status": [(line_id, 1)]})
        return action

agent = ReconnectAgent(env.action_space)
obs = env.reset()
reward = env.reward_range[0]
done = False
total_reward = 0

for step in range(env.chronics_handler.max_timestep()):
    action = agent.act(obs, reward, done)
    obs, reward, done, info = env.step(action)
    total_reward += reward
    if done:
        break

print(f"Survived {step} steps, total reward: {total_reward:.2f}")

The L2RPN (Learning to Run a Power Network) competition uses Grid2Op and has driven significant advances in ML-based grid control. Winning agents combine imitation learning from expert operators with RL fine-tuning.

Capacity expansion planning with PyPSA

For long-term planning (deciding which generators, storage, and lines to build):

import pypsa

network = pypsa.Network()
network.set_snapshots(pd.date_range("2025-01-01", periods=8760, freq="1h"))

# Add bus
network.add("Bus", "main_bus")

# Add existing gas plant
network.add("Generator", "gas", bus="main_bus",
            p_nom=100, marginal_cost=50,  # $/MWh
            carrier="gas")

# Add candidate solar (optimizable capacity)
network.add("Generator", "solar", bus="main_bus",
            p_nom_extendable=True, p_nom_max=500,
            capital_cost=1200000,  # $/MW annualized
            marginal_cost=0,
            p_max_pu=solar_capacity_factor_series,
            carrier="solar")

# Add candidate battery
network.add("StorageUnit", "battery", bus="main_bus",
            p_nom_extendable=True, p_nom_max=200,
            capital_cost=300000,  # $/MW annualized
            max_hours=4,  # 4-hour duration
            efficiency_store=0.95, efficiency_dispatch=0.95,
            carrier="battery")

# Add load
network.add("Load", "demand", bus="main_bus",
            p_set=load_profile_series)

# Solve capacity expansion
network.optimize(solver_name="glpk")

print(f"Optimal solar capacity: {network.generators.p_nom_opt['solar']:.0f} MW")
print(f"Optimal battery capacity: {network.storage_units.p_nom_opt['battery']:.0f} MW")
print(f"Total system cost: ${network.objective:.0f}")

Tradeoffs

DecisionOption AOption B
Solverpandapower (detailed distribution)PyPSA (simplified, large-scale optimization)
Time resolutionHourly (standard, faster)15-minute (captures fast transients, 4× more data)
Control approachRule-based (predictable, suboptimal)RL-based (adaptive, harder to verify safety)
Network detailFull 3-phase unbalancedSingle-phase equivalent (faster, less accurate for LV networks)
Planning horizon1 year operational (detailed)20+ year capacity expansion (strategic)

One thing to remember: Smart grid simulation is where physics meets optimization — the power flow equations define what’s physically possible, and optimization algorithms find the best operating point within those constraints. Python bridges both worlds.

pythonsmart-gridenergysimulation

See Also

Related Topics

ab-testing How tech companies run thousands of experiments at once to improve their products — the scientific method applied to everything from button colors to recommendation algorithms. 3 levels → activation-functions Why neural networks need these tiny mathematical functions — and how ReLU's simplicity accidentally made deep learning possible. 3 levels → ai-agents-architecture How AI systems go from answering questions to actually doing things — the design patterns that turn language models into autonomous agents that browse, code, and plan. 3 levels →