Analytical Uncertainty Propagation for Multi-Period Stochastic Optimal Power Flow

How can we combine stochastic DC optimal power flow with Gaussian processes?

This repository provides both the finished paper (ArXiv) and a Julia package to tackle this question.

Developing with Julia

Getting started with Julia: getting started, workflow tips, or notable differences from other languages

To develop a Julia package the package Revise (using Revise) is helpful, as it allows changes in a file to be applied directly, so that you don't have to restart Julia. Use by using Revise.

Many questions are answered in Julias discourse forum.

Installation

The code should be self-contained and run out of the box.

First, install Julia (we used Julia 1.3).

Packages

Required packages: PyPlot, JuMP, MosekTools, PowerModels, MAT, LinearAlgebra, StatsFuns. To install packages

• start Julia
• enter ] to switch the the package manager. Your console should read (v1.3) pkg> (depends on your Julia version)
• enter add PyPlot JuMP MosekTools PowerModels MAT StatsFuns.

Additionally, the code requires two self-written Julia packages that need to be registered locally. Note that this very repository is a package, namely the Julia Package DCsOPF.

• clone the repo for PowerModelsParsing to <your-path>/PowerModelsParsing
• clone this repo to <your-path>/DCsOPF
• start Julia
• enter ] to switch the the package manager. Your console should read (v1.3) pkg> (depends on your Julia version)
• enter dev <your-path>/PowerModelsParsing
• enter dev <your-path>/DCsOPF.

For more information on adding unregistered packages, see the Julia package manager documentation.

Mosek license

The code needs Mosek to solve the optimization problem. You can obtain a license here. The Julia package MosekTools will complain if it can't find the license file, and it should tell you where to put it by default. In case of any issues, consult the MosekTools documentation

Run & modify the code

Run the code

To run the code simply execute the main.jl file from top to bottom. At the top there are parameters for setup of plots, then there are different cases, and then the main pipeline.

Modify the code

We can modify the code by adding a new case (casefiles/case.m) and by adding new wind data (data/CovMatrix.mat).

To add a new case (network) the MatPower casefile case.m with the network topology and parameters is needed, as well as the PTDF (power transfer distribution factor) file. The wind data (covariance matrix) can be applied to any case.

Description	Example	Note
Matpower case file	casefiles/case5.m	Has to be of a specific caseformat e.g. IEEE case5
Power transfer distribution factor matrix	casefiles/case5_PTDF.mat	The PTDF can be extracted out of case.m with MATLAB using mpc = loadcase('<case>.m'); ptdf = makePTDF(mpc); save('<case>_PTDF.mat'), 'ptdf');.
GP wind data (covariance matrix)	data/CovMatrix.mat	The GP (Gaussian process) is created with GPR (Gaussian process regression), e.g. with the GPflow package.

Gaussian processes (covariance matrices)

The wind forecast is given in form of a covariance matrix, e.g. CovMatrix_artificial.mat, that contains the mean and variance of the forecast. If we forecast over n time steps, then the files are:

• mu_post: mean (size nx1)
• Lpost: variance (size nxn, lower triangular matrix) Letting Σ = Lpost * Lpost', then (mu_post, Σ) is a Gaussian process over n time steps.

Add case in main.jl

To test the code with your own network add the parameters for another case with the following code:

params["<case>"] = CaseParams("<case number>", # case number, e.g. "5"
unc_set, # wind buses (uncertainties), e.g. [[1,2],[1,2,3]]
stor_set, # storage busses, e.g. [[4],[4,5]]
true, # executed local optimization
true, # execute global optimization
false, # solve deterministic problem (covariance matricies are zero)
false, # use artificial wind power curves
1., # multiply load with factor
1., # multiply wind injection with factor
10., # multiply uncertainty (variance) of wind with factor
10, # storage upper bound
0.8) # amplitude of load modelled as sine curve

where each load is a dictionary with two fields :μ and :Σ. This is where you can insert your GP information. Just make sure that the numerical range is about the same. Also note that loads are negative injections, hence the minus sign.

GPR

For the Gaussian process regression we use GPflow. These files are not currently for display.

Data

The wind data is taken from Zenodo containing hourly data from 2014 to 2021.

Questions?

Write to us!