## Optimal Power Flow
_**[Power Systems Optimization](https://github.com/east-winds/power-systems-optimization)**_

_by Michael R. Davidson, Jesse D. Jenkins, and Sambuddha Chakrabarti_

This notebook consists an introductory glimpse of and a few hands-on activities and demostrations of the Optimal Power Flow (OPF) problem—which minimizes the short-run production costs of meeting electricity demand from a given set of generators subject to various technical and flow limit constraints.

We will talk about a single-time period, simple generator, and line flow limit constraints (while modeling the network flows as dictated by the laws of physics). This is adds a layer of complexity and sophistication on top of the Economic Dispatch (ED) problem.

Since we will only discuss single time-period version of the problem, we will not be considering inter-temporal constraints, like ramp-rate limits. However, this model can easily be extended to allow for such constraints.

We will start off with some simple systems, whose solutions can be worked out manually without resorting to any mathematical optimization model and software. But, eventually we will be solving larger system, thereby emphasizing the importance of such software and mathematical models.

## Introduction to OPF

Optimal Power Flow (OPF) is a power system optimal scheduling problem which fully captures the physics of electricity flows, which adds a leyr of complexity, as well gives a more realistic version of the Economic Dispatch (ED) problem. It usually attempts to capture the entire network topology by representing the interconnections between the different nodes through transmission lines and also representing the electrical parameters, like the resistance, series reactance, shunt admittance etc. of the lines. however, the full-blown "AC" OPF turns out to be an extremely hard problem to solve (usually NP-hard). Hence, system operators and power marketers usually go about solving a linearized version of it, called the DC-OPF. The DC-OPF approximation works satisfactorily for bulk power transmission networks as long as such networks are not operated at the brink of instability or, under very heavily heavily loaded conditions.

## Single-time period, simple generator constraints
We will first examine the case where we are optimizing dispatch for a single snapshot in time, with only very simple constraints on the generators. $x^2$


$$
\begin{align}
		\mathbf{Objective\;Function:}\min_{P_g}\sum_{g\in{G}}C_{g}(P_{g})\longleftarrow\mathbf{power\;generation\; cost}\\
		\mathbf{Subject\;to:\:}{\underline{P}_{g}}\leqslant{P_{g}}\leqslant{{\overline{P}_{g}}},\;\forall{g\in{G}}\longleftarrow\mathbf{MW\; generation\; limits}\\
		P_{g(i)}-P_{d(i)}\longleftarrow\mathbf{real\; power\; injection}\notag\\=\sum_{j\in J(i)}B_{ij}(\theta_j-\theta_i),\;\forall{{i}\in\mathcal{N}}\\
		|P_{ij}|\leqslant{\overline{P}_{ij}},\;\forall{ij}\in{T}\longleftarrow\mathbf{MW\; line\; limit}\\
\end{align}
$$
 
The **decision variable** in the above problem is:

- $P_{g}$, the generation (in MW) produced by each generator, $g$
- $\theta_i$, $\theta_j$ the voltage phase angle of each bus/node, $i,j$

The **parameters** are:

- ${\underline{P}_{g}}$, the minimum operating bounds for the generator (based on engineering or natural resource constraints)
- ${\overline{P}_{g}}$, the maximum operating bounds for the generator (based on engineering or natural resource constraints)
- $P_{d(i)}$, the demand (in MW) at node $i$
- ${\overline{P}_{ij}}$, the line-flow limit for line connecting buses $i$ and $j$
- $B_{ij}$, susceptance for line connecting buses $i$ and $j$

just like the ED problem, here also, we can safely ignore fixed costs for the purposes of finding optimal dispatch.

With that, let's implement OPF.

# 1. Load packages¶


```julia
# New packages introduced in this tutorial (uncomment to download the first time)
import Pkg; Pkg.add("PlotlyBase")
using JuMP, GLPK
using Plots; plotly();
using VegaLite  # to make some nice plots
using DataFrames, CSV, PrettyTables
ENV["COLUMNS"]=120; # Set so all columns of DataFrames and Matrices are displayed
```

    [32m[1m  Updating[22m[39m registry at `~/.julia/registries/General`
    [32m[1m  Updating[22m[39m git-repo `https://github.com/JuliaRegistries/General.git`
    [2K[?25h[1mFetching:[22m[39m [========================================>]  100.0 %.0 %13.1 %>                             ]  26.1 %[36m[1mFetching:[22m[39m [================>                        ]  39.2 %                   ]  52.2 % [===========================>             ]  65.2 %78.2 %[32m[1m Resolving[22m[39m package versions...
    [32m[1m Installed[22m[39m LaTeXStrings ──────── v1.2.0
    [32m[1m Installed[22m[39m PlotlyBase ────────── v0.4.1
    [32m[1m Installed[22m[39m DocStringExtensions ─ v0.8.3
    [32m[1m  Updating[22m[39m `~/.julia/environments/v1.3/Project.toml`
     [90m [a03496cd][39m[92m + PlotlyBase v0.4.1[39m
    [32m[1m  Updating[22m[39m `~/.julia/environments/v1.3/Manifest.toml`
     [90m [ffbed154][39m[92m + DocStringExtensions v0.8.3[39m
     [90m [b964fa9f][39m[92m + LaTeXStrings v1.2.0[39m
     [90m [a03496cd][39m[92m + PlotlyBase v0.4.1[39m


    ┌ Info: Precompiling PlotlyBase [a03496cd-edff-5a9b-9e67-9cda94a718b5]
    └ @ Base loading.jl:1273


### 2. Load and format data

We will use data for IEEE 118 bus test case and two other test cases for a 3 bus and a 2 bus system:

- generator cost curve, power limit data, and connection-node
- load demand data with MW demand and connection node
- transmission line data with resistance, reactance, line MW capacity, from, and to nodes


```julia
datadir = joinpath("OPF_data") 
# Note: joinpath is a good way to create path reference that is agnostic
# to what file system you are using (e.g. whether directories are denoted 
# with a forward or backwards slash).
gen_info = CSV.read(joinpath(datadir,"Gen118.csv"), DataFrame);
line_info = CSV.read(joinpath(datadir,"Tran118.csv"), DataFrame);
loads = CSV.read(joinpath(datadir,"Load118.csv"), DataFrame);

# Rename all columns to lowercase (by convention)
for f in [gen_info, line_info, loads]
    rename!(f,lowercase.(names(f)))
end
```


```julia
#=
Function to solve Optimal Power Flow (OPF) problem (single-time period)
Inputs:
    gen_info -- dataframe with generator info
    line_info -- dataframe with transmission lines info
    loads  -- dataframe with load info
Note: it is always a good idea to include a comment blog describing your
function's inputs clearly!
=#
function OPF_single(gen_df, line_info, loads)
    OPF = Model(GLPK.Optimizer) # You could use Clp as well, with Clp.Optimizer
    
    # Define sets based on data
      # A set of all variable generators
    G_var = gen_df[gen_df[!,:is_variable] .== 1,:r_id] 
      # A set of all non-variable generators
    G_nonvar = gen_df[gen_df[!,:is_variable] .== 0,:r_id]
      # Set of all generators
    G = gen_df.r_id
    # Extract some parameters given the input data
      # Generator capacity factor time series for variable generators
    gen_var_cf = innerjoin(gen_variable, 
                    gen_df[gen_df.is_variable .== 1 , 
                        [:r_id, :gen_full, :existing_cap_mw]], 
                    on = :gen_full)
        
    # Decision variables   
    @variables(ED, begin
        GEN[G]  >= 0     # generation
        # Note: we assume Pmin = 0 for all resources for simplicty here
    end)
                
    # Objective function
    @objective(ED, Min, 
        sum( (gen_df[i,:heat_rate_mmbtu_per_mwh] * gen_df[i,:fuel_cost] +
            gen_df[i,:var_om_cost_per_mwh]) * GEN[i] 
                        for i in G_nonvar) + 
        sum(gen_df[i,:var_om_cost_per_mwh] * GEN[i] 
                        for i in G_var)
    )

    # Demand constraint
    @constraint(ED, cDemand, 
        sum(GEN[i] for i in G) == loads[1,:demand])

    # Capacity constraint (non-variable generation)
    for i in G_nonvar
        @constraint(ED, GEN[i] <= gen_df[i,:existing_cap_mw])
    end

    # Variable generation capacity constraint
    for i in 1:nrow(gen_var_cf)
        @constraint(ED, GEN[gen_var_cf[i,:r_id] ] <= 
                        gen_var_cf[i,:cf] *
                        gen_var_cf[i,:existing_cap_mw])
    end

    # Solve statement (! indicates runs in place)
    optimize!(ED)

    # Dataframe of optimal decision variables
    solution = DataFrame(
        r_id = gen_df.r_id,
        resource = gen_df.resource,
        gen = value.(GEN).data
        )

    # Return the solution and objective as named tuple
    return (
        solution = solution, 
        cost = objective_value(ED),
    )
end
```




    OPF_single (generic function with 1 method)




```julia

```

a14411ae2d618522ce56b738a289f41f944a5d5f
ipynb
Jupyter Notebook
Notebooks/06-OPF-problem_other.ipynb
sambuddhac/power-systems-optimization
f65ec4b718807703452cf6723105926dac73c649
CC-BY-4.0 MIT
Notebooks/06-OPF-problem_other.ipynb
sambuddhac/power-systems-optimization
f65ec4b718807703452cf6723105926dac73c649
CC-BY-4.0 MIT
Notebooks/06-OPF-problem_other.ipynb
sambuddhac/power-systems-optimization
f65ec4b718807703452cf6723105926dac73c649
CC-BY-4.0 MIT
Qwen/Qwen-72B
1. YES
2. YES