Instanton demo

The function instantonAnalysis handles several varieties of the instanton problem:

instantonAnalysis(psData; method="lagrange",
                          constraintType="activeFlow",
                          genResponse="droop")

The three keyword arguments determine which instanton analysis method to call.

Let's use Ipopt to solve the instanton problem with active power flows and proportional conventional generator response. First, make sure all necessary packages are available:

Now make sure the julia-code directory is beneath the current directory:

Now we include instanton Julia code:

Load power system data from a MATLAB .mat file into an instance of custom powerSystemData type, then perform instanton analysis using JuMP with Ipopt:

To display results, import SVG display capability for the notebook, then call plotInstantonResults:

How it works: the JuMP model

Loading data from the .mat file into a custom Julia type is straightforward. The type has a field for each relevant piece of power system data:

Similarly, I define a type to store analysis results:

Having stored all system data in a single object, I pass it to the following function, which builds a series of models, passes them to Ipopt, and returns solutions. Here is an abbreviated version of the code. In a moment we will break it down a few lines at a time.

Let's concentrate on the solver model composed with JuMP macros.

1. Create a new model instance:

   m = Model(solver = IpoptSolver())
   

2. Define some variables:

   N = length(Rp)
   @defVar(m, ρ[1:N] >= 0) # Add decision variables
   @defVar(m, θ[1:N]) # Add bus angles
   @defVar(m, α) # mismatch
   

   Compare with:

   • $\rho_i$ is renewable generation at bus $i$ in per unit (non-negative).
   • $\theta_k$ is the phase angle at bus $k$.
   • $\alpha$ is the participation coefficient (mismatch), defined as $$\alpha:= \sum D - \sum \rho - \sum G_0 ~.$$

3. Set the objective:

   setObjective(m, :Min, 0.5*sum([(ρ[i] - Rp[i])^2 
                       for i in find(Rp)]))
   

   Compare with: $$\min \frac{1}{2} \left( \rho - \rho_0 \right)^\top \left( \rho - \rho_0 \right)$$

4. Power balance constraints:

   @defConstrRef PowerBalance[1:N]
   for i in setdiff(1:N,find(bustype.==3))
    PowerBalance[i] = @addConstraint(m, 
        sum([Y[i,j]*θ[j] for j in 1:N]) == 
        Gp[i] + k[i]*α + ρ[i] - Dp[i])
   end
   

   Compare with: $$\begin{align} \sum_k( Y_{ik} \theta_k) &= G_{i} + \rho_i - D_{i} \quad \forall i \in N \\ G &= G_{0} + k\alpha \end{align}$$

5. Constrain renewable generation to 0 for non-wind nodes:

   @addConstraint(m, NonWind, 
    sum([ρ[i] for i in setdiff(collect(1:N),find(Rp))]) == 0)
   

6. Define slack bus:

   @addConstraint(m, Slack, 
    θ[find(bustype.==3)[1]] == 0)
   

   Compare with: $$\theta_i = 0 \quad \text{ where bus $i$ is the angle reference} $$

7. Add "congestion constraint":

   @addConstraint(m, Congestion, 
    θ[f[idx]] - θ[t[idx]] == sense*x[idx]*Plim[idx])
   

   (Note that sense is $\pm 1$.) Compare with: $$ \theta_i - \theta_k = x_{ik} P_{lim,ik} \quad \text{for each } (i,k) \in G $$

8. Conventional generator droop response:

   @addConstraint(m, Participation,
    α == (sum(Dp) - sum([ρ[i] for i in find(Rp)])) - sum(Gp))
   

   Compare with:

   • $\alpha$ is the participation coefficient (mismatch), defined as

For each line (and direction) in the network, a model is generated and passed to Ipopt. The last line of the solver_activeFlow_droop function returns solution information:

return m, ρ, θ, α, PowerBalance, Slack, Congestion, Participation

How it works: displaying results

The plotInstantonResults() function generates the graph we saw:

function plotInstantonResults(  eventIdx,
                                psData,
                                results;
                                plotType="graph",
                                plotName="defaultPlotName",
                                constraintType="activeFlow")

Note the presence of positional and keyword arguments. If we want a visualization other than the graph, we can choose another value for plotType and add code to the function to respond to it. The function is general enough to handle all plotting.

The graph output in particular uses a function I wrote called writeDot(), which is about 70 lines of code and generates a .dot file. Graphviz uses this file to render the network as a directed graph.

What we've covered so far

Here we ran Julia code to perform instanton analysis on the IEEE RTS-96 network. We stepped through the process of translating a mathematical programming model into Julia JuMP macros. We also plotted the results using Graphviz and discussed how this was implemented in Julia.

Where to next?

In the last notebook we will explore an extension to the instanton analysis presented here. Along the way we will see the usefulness of Julia and the IJulia notebook in exploration, prototyping, and presentation.