Library

FlexibilityAnalysis

A Julia package that implements a computational framework for quantifying and analyzing system flexibility.

The basic functionality emulates typical JuMP models to facilitate general computation of the flexibility index problem. Thus, basic functionality uses FlexibilityModel, @randomvariable, @recoursevariable, @variable, @constraint, setuncertaintyset, and solve.

Methods/Macros

• FlexibilityModel
• @randomvariable
• @recoursevariable
• getflexibilitydata
• setcovariance
• getcovariance
• getmean
• setmean
• setuncertaintyset
• ismeanfeasible
• findcenteredmean
• getflexibilityindex
• getsolutiontime
• getconfidencelevel
• getactiveconstraints
• rankinequalities
• findstochasticflexibility

FlexibilityModel(; [solver = JuMP.UnsetSolver()])

Return a flexibility model object which extends a JuMP model object to contain FlexibilityData and implement a custom solvehook. An appropriate solver should be specified in order solve the flexibility index problem. A solver capable of handling MIQCPs is required for ellipsoidal and 2-norm uncertainty sets otherwise a MILP solver can be used. This model is solved with solve, see solvehook for documention on the accepted keyword arguments.

Arguments

• solver = JuMP.UnsetSolver() The solver, should use an MIQCP, MINLP, or MILP solver as appropriate.

julia> m = FlexibilityModel(solver = GurobiSolver())
Feasibility problem with:
 * 0 linear constraints
 * 0 variables
Solver is Gurobi

@randomvariable(m, x, mean)

Defines a random variable using RandomVariable(m::Model, mean::Number, name::AbstractString) and requires that a mean for the variable be provided. This can later be overwritten with setmean.

Arguments

• m::Model The flexibility model.
• x::Symbol The variable name.
• mean::Number The variable mean

julia> @randomvariable(m2, w, mean = 42)

julia> @randomvariable(m2, ws[i = 1:4], mean = 42)
4-element Array{FlexibilityAnalysis.RandomVariable,1}:
 FlexibilityAnalysis.RandomVariable(Feasibility problem with:
 * 9 linear constraints
 * 89 variables

 julia> @randomvariable(m2, ws[i = 1:4], mean = [1; 2; 3; 4][i])
 4-element Array{FlexibilityAnalysis.RandomVariable,1}:
  FlexibilityAnalysis.RandomVariable(Feasibility problem with:
  * 9 linear constraints
  * 96 variables

@recoursevariable(m, x)

Defines a recourse variable using RecourseVariable(m::Model, name::AbstractString).

Arguments

• m::Model The flexibility model.
• x::Symbol The variable name.

julia> @recoursevariable(m2, d)

julia> @recoursevariable(m2, ds[1:4])
4-element Array{FlexibilityAnalysis.RecourseVariable,1}:
 FlexibilityAnalysis.RecourseVariable(Feasibility problem with:
 * 9 linear constraints
 * 101 variables

getflexibilitydata(m::Model)

Return the FlexibilityData corresponding the flexibility model m. An error is thrown if m is a regular JuMP model.

Arguments

• m::Model The flexibility model.

julia> getflexibilitydata(m)
FlexibilityAnalysis.FlexibilityData(JuMP.AbstractConstraint[], 0, Number[], AbstractString[], Int64[], 0, AbstractString[], Int64[], FlexibilityAnalysis.EllipsoidalSet(:Ellipsoid, false), Array{Number}(undef,0,0), nothing, Int64[], nothing)

setcovariance(m::Model, covariance::Matrix)

Specify the covariance matrix covariance to be stored in the flexibility model m. This method verifies that the matrix is symmetric positive semi-definite and writes it to FlexibilityData.covariance.

Arguments

• m::Model The flexibility model.
• covariance::Matrix The covariance matrix.

julia> setcovariance(m, [2 1; 1 2])
2×2 Array{Int64,2}:
 2  1
 1  2

getcovariance(m::Model)

Return the current covariance matrix in the flexibility model as stored in FlexibilityData.covariance.

Arguments

• m::Model The flexibility model.

julia> getcovariance(m)
4×4 Array{Number,2}:
 11.11   0.0    0.0    0.0
  0.0   11.11   0.0    0.0
  0.0    0.0   11.11   0.0
  0.0    0.0    0.0   11.11

setmean(m::Model, mean::Vector)

Specify the mean corresponding to FlexibilityData.RVmeans stored in the flexibility model m. This method verifies that the length of the input mean matches the length of FlexibilityData.RVmeans before overwriting the current mean.

Arguments

• m::Model The flexibility model.
• mean::Vector The means of the random variables.

setmean(m, [2.3; 5])
2-element Array{Float64,1}:
 2.3
 5.0

getmean(m::Model)

Return the current mean vector in the flexibility model as stored in FlexibilityData.RVmeans.

Arguments

• m::Model The flexibility model.

julia> getmean(m)
4-element Array{Number,1}:
 620
 388
 583
 313

getmean(variable::RandomVariable)

Return the mean corresponding to a particular RandomVariable. Currently this only accepts a single random variable and vector variables are not accepted directly.

Arguments

• variable::RandomVariable The random variable name, must be a single variable.

julia> getmean(T[1])
620

setuncertaintyset(m::Model, uncertainty_set::Symbol, [attribute = nothing; only_positive::Bool = false])

Specify the type of uncertainty set to be used and stored in FlexibilityData.uncertainty_set and provide the necessary attribute. An ellipsoidal uncertainty set can be specified with the :Ellipsoid symbol and the corresponding covariance matrix will need to input via attribute if it has not already been set with setcovariance. A hyperbox uncertainty set is specified with the :Hyperbox symbol and corresponding negative/positive deviation vectors need to be inputed via attribute as a vector of vectors of the form [[neg_dev]; [pos_dev]]. Finally, a p-norm uncertainty set can be specified with the :PNorm symbol and providing the corresponding p value via attribute where p can equal 1, 2, or Inf.

Arguments

• m::Model The flexibility model.
• uncertainty_set::Symbol The uncertainty set name.
• attribute = nothing The necessary atribute for the specified uncertainty set.

Keyword Arguments

• only_positive::Bool = false Indicate if the uncertainty set should be intersected with the set of all positive values.

julia> setuncertaintyset(m, :Ellipsoid)

julia> setuncertaintyset(m, :Ellipsoid, eye(4)* 11.11)

julia> setuncertaintyset(m, :Hyperbox, [[ones(4)]; [ones(4)]])
FlexibilityAnalysis.HyperboxSet(:Hyperbox, Number[1.0, 1.0, 1.0, 1.0], Number[1.0, 1.0, 1.0, 1.0], false)

julia> setuncertaintyset(m, :PNorm, 1)
FlexibilityAnalysis.PNormSet(:PNorm, 1, false)

julia> setuncertaintyset(m, :PNorm, 2, only_positive = true)
FlexibilityAnalysis.PNormSet(:PNorm, 2, true)

julia> setuncertaintyset(m, :PNorm, Inf)
FlexibilityAnalysis.PNormSet(:PNorm, Inf, false)

ismeanfeasible(m::Model; [toler::Number = 1e-5, solver = Clp.ClpSolver()])

Returns a Bool indicating if the mean stored in FlexibilityData.RVmeans is feasible, meaning that it lies inside the feasible region. This check is done using the so-called feasibility function.

Arguments

• m::Model The flexibility model.

Keyword Arguments

• toler::Number = 1e-5 The numerical tolerance for checking the feasibility.
• solver = Clp.ClpSolver() The solver, any LP or NLP solver shoudl work.

julia> ismeanfeasible(m)
true

julia> ismeanfeasible(m, solver = GurobiSolver(OutputFlag = 0))
Academic license - for non-commercial use only
true

findcenteredmean(m::Model; [center::Symbol = :feasible, solver = Ipopt.IpoptSolver(print_level = 0), toler::Number = 1e-5, update_mean::Bool = false, only_positive::Bool = false])

Returns a center point based on the analytic or feasible center. The result can overwrite the mean stored in FlexibilityData.RVmeans if desired. This is a wrapper function for ComputeCenter.

Arguments

• m::Model The flexibility model.

Keyword Arguments

• center::Symbol = :feasible Indicates the type of center, accepted arguments are :feasible and :analytic.
• solver = Ipopt.IpoptSolver(print_level = 0) The solver which must be an NLP solver for the analytic center.
• toler::Number = 1e-5 The tolerance to check solution validity.
• update_mean::Bool = false Indicates if the computed center should overwrite FlexibilityData.RVmeans.
• only_positive::Bool = false Indicates if the center need by strictly positive.

julia> findcenteredmean(m, only_positive = true)
4-element Array{Float64,1}:
 1684.74
   79.0718
  195.073
    0.0

julia> findcenteredmean(m, center = :analytic, update_mean = true)
4-element Array{Float64,1}:
  898.125
 -507.214
  594.544
  317.23

getflexibilityindex(m::Model)

Return the current flexibility index in the flexibility model as stored in FlexibilityData.flexibility_index.

Arguments

• m::Model The flexibility model.

julia> getflexibilityindex(m)
3.5993764186390327

getconfidencelevel(m::Model)

Return the confidence level provided that the flexibility model m has been solved with an ellipsoidal uncertainty set. This is equivalent to the quantile associated with the optimized uncertainty set. Note that this assumes a multivariate Gaussian distribution with mean FlexibilityData.RVmeans and covariance FlexibilityData.covariance.

Arguments

• m::Model The flexibility model.

julia> getconfidencelevel(m)
0.5370703369769008

getsolutiontime(m::Model)

Return the solution time to compute the flexibility index as stored in the flexibility model as stored in FlexibilityData.solution_time.

Arguments

• m::Model The flexibility model.

julia> getsolutiontime(m)
0.00199127197265625

getactiveconstraints(m::Model)

Return the current vector of active constraint indexes in the flexibility model as stored in FlexibilityData.active_constraints.

Arguments

• m::Model The flexibility model.

julia> getactiveconstraints(m)
2-element Array{Int64,1}:
 3
 6

JuMP.getvalue(v::FlexibilityVariable)

Return the value of the a flexibility variable this is an extension of JuMP.getvalue.

Arguments

• v::FlexibilityVariable The flexibility variable.

julia> getvalue(T)
4-element Array{Float64,1}:
 620.0
 388.0
 581.0
 319.0

rankinequalities(m::Model; max_ranks::Int = 5, suppress_warnings::Bool = true, U::Int = 10000, diag::Bool = false, active_constr::Bool = false, real_recourse_dim::Int = -1, conic_δ::Bool = false)

Returns ranking data in the form Vector{Dict} where each dictionary corresponds to a particular rank level and contains the optimal flexibility_index, active constraint indexes, and flexibility model. The function will iteratively solve copies of the flexibility model via solve where the prior active constraints are turned off in order rank the most limiting constraints. The user can specify the maximum number of rank levels and the flexibility index problem will be repeatedly solved until that maximum is acheived or the problem becomes unbounded, which occurs first.

Arguments

• m::Model The flexibility model.

Keyword Arguments

• max_ranks::Int = 5 The maximum number of rank levels. 2
• suppress_warnings::Bool = false Indicates if solver warnings should be suppressed.
• U::Union{Int, Float64} = 10000 The slack variable upper bound.
• diag::Bool = false Indicates whether or not to diagnonalize ellipsoidal uncertainty set (this is only active when an ellipsoidal set is used).
• active_constr::Bool = false Indicates if the optional active constraint should be used which enforces how many inequalities are active at the solution, this must be set to true for systems without control variables and/or contain state variables.
• real_recourse_dim::Int = -1 The actual number of recourse variables in case state variables are included as recourse variables. This is mandatory if active_constr = true and no state variables are provided.
• conic_δ::Bool = false This should be set to true if a conic solver is used such as Pajarito.jl.

julia> rankinequalities(m, active_constr = true)
2-element Array{Dict,1}:
 Dict{String,Any}(Pair{String,Any}("flexibility_index", 3.59938),Pair{String,Any}("model", Feasibility problem with:
 * 6 linear constraints
 * 6 variables
Solver is Pavito),Pair{String,Any}("active_constraints", [3, 6]))
 Dict{String,Any}(Pair{String,Any}("flexibility_index", 9.58983),Pair{String,Any}("model", Feasibility problem with:
 * 6 linear constraints
 * 6 variables
Solver is Pavito),Pair{String,Any}("active_constraints", [1, 5]))

findstochasticflexibility(m::Model; num_pts::Int = 10000, toler::Number = 1e-5, solver = Clp.ClpSolver(), only_positive::Bool = false, use_vulnerability_model::Bool = false, use_flexibility_index::Bool = false, seed::Int = -1)

Returns the estimated stochastic flexibility index that is evaluated via Monte Carlo sampling. At default this estimation is carried out by evaluating the feasibility of each Monte Carlo sample. The samples are generating from a multivariate Gaussian distribution with mean FlexibilityData.RVmeans and covariance FlexibilityData.covariance. The vulnerability model also tests the feasibility of the samples, but does so in one large optimization problem instead of evaluating each sample indiviually. The optimized flexibility index can also be used to reduce the number of samples that need to be evaluated.

Arguments

• m::Model The flexibility model.

Keyword Arguments

• num_pts::Int = 10000 Number of Monte Carlo samples.
• toler::Number = 1e-5 The feasibility check tolerance.
• solver = Clp.ClpSolver() The solver, any LP or NLP solver should work.
• only_positive::Bool = false Indicates if only positive samples should be used.
• use_vulnerability_model::Bool = false Indicates if the vulnerability model should be used.
• use_flexibility_index::Bool = false Indicates if the optimal flexibility index should be used.
• seed::Int = -1 Random seed for sample collection, any negative value will turn off the random seed.

julia> findstochasticflexibility(m)
0.9687

julia> findstochasticflexibility(m, use_vulnerability_model = true)
0.9705

julia> findstochasticflexibility(m, num_pts = 5000, use_flexibility_index = true)
0.973

JuMP.linearindex(v::FlexibilityVariable)

Return the index of the a flexibility variable this is an extension of JuMP.linearindex.

Arguments

• v::FlexibilityVariable The flexibility variable, must be a single variable.

julia> linearindex(Qc)
1

FlexibilityVariable <: JuMP.AbstractJuMPScalar

An abstract type to define new variable types.

RandomVariable <: FlexibilityVariable

A DataType for random variables.

Fields

• m::Model Flexibility model.
• idx::Int Index of variable in model.

RecourseVariable <: FlexibilityVariable

A DataType for recourse variables.

Fields

• m::Model Flexibility model.
• idx::Int Index of variable in model.

FlexibilityExpr <: JuMP.GenericAffExpr

A GenericAffExpr that contains random and/or recourse variables.

FlexibilityConstraint <: JuMP.AbstractConstraint

A constraint that contains random and/or recourse variables.

Fields

• flex_expr::FlexibilityExpr Constraint expression.
• sense::Symbol The cosntraint sense symbol :(<=) or :(>=) or :(==).

AbstractUncertaintySet

An abstract type to define new uncertainty set types.

EllipsoidalSet <: AbstractUncertaintySet

An ellipsoidal uncertainty set that will use the covariance matrix stored in FlexibilityData.

Fields

• name::Symbol Name of the set which will be :Ellipsoid.
• only_positive::Bool An option to indicate if the set should be intersected with the set all positive numbers.

HyperboxSet <: AbstractUncertaintySet

A hyperbox uncertainty set whose nomimal dimensions are determined by neg_dev and pos_dev.

Fields

• name::Symbol Name of the set which will be :Hyperbox.
• neg_dev::Vector{Number} A vector of the expected negative deviation of the random variables.
• pos_dev::Vector{Number} A vector of the expected positive deviation of the random variables.
• only_positive::Bool An option to indicate if the set should be intersected with the set all positive numbers.

PNormSet <: AbstractUncertaintySet

A p-norm based uncertainty set based on a bounded p-norm.

Fields

• name::Symbol Name of the set which will be :PNorm.
• p::Number The value of p which can be 1, 2, or Inf.
• only_positive::Bool An option to indicate if the set should be intersected with the set all positive numbers.

FlexibilityData

A DataType for storing the data necessary to manage the bookkeeping of the flexibility variables (RandomVariable and RecourseVariable), the uncertainty set, and solution results.

Fields

• flexibility_constraints::Vector{JuMP.AbstractConstraint} Constraints that involve flexibility variables.
• numRVs::Int The number of RandomVariable that have been added to the model.
• RVmeans::Vector{Number} The means corresponding to each RandomVariable.
• RVnames::Vector{AbstractString} The symbolic name of each RandomVariable.
• RVcols::Vector{Int} The index of each RandomVariable.
• num_recourse_vars::Int The number of RecourseVariable that have been added to the model.
• recourse_names::Vector{AbstractString} The symbolic name of each RecourseVariable.
• recourse_cols::Vector{Int} The index of each RecourseVariable.
• uncertainty_set::AbstractUncertaintySet The uncertainty set DataType with all of the set specfic attributes.
• covariance::Matrix{Number} The covariance matrix.
• flexibility_index::Union{Nothing, Number} The flexibility index result obtained from solving the flexibility model.
• active_constraints::Vector{Int} The indexes of the active inequality constraints at the solution of the flexibility model.
• 'solution_time::Number' The solution time in seconds.

getuncertaintyset(m::Model)

Return the current uncertainty set datatype in the flexibility modelas stored in FlexibilityData.uncertaintyset.

Arguments

• m::Model The flexibility model.

solvehook(m::Model; [suppress_warnings::Bool = false, U::Number = 10000, diag::Bool = false, active_constr::Bool = false, real_recourse_dim::Int = -1, conic_δ::Bool = false, inactives::Vector = []])

Returns the solution status to solving the flexibility model m. This solvehook what JuMP.solve(::Model) for flexibility models. This solves the flexibility index problem using the variables and constraints specified in m.

Arguments

• m::Model The flexibility model.

Keyword Arguments

• suppress_warnings::Bool = false Indicates if solver warnings should be suppressed.
• U::Number = 10000 The slack variable upper bound.
• diag::Bool = false Indicates if the ellipsoidal uncertainty set is diagonalized (this is only active when an ellipsoidal set is used).
• active_constr::Bool = false Indicates if the optional active constraint should be used which enforces how many inequalities are active at the solution, this must be set to true for systems without control variables and/or contain state variables.
• real_recourse_dim::Int = -1 The actual number of recourse variables in case state variables are included as recourse variables. This is mandatory if active_constr = true and no state variables are provided.
• conic_δ::Bool = false This should be set to true if a conic solver is used such as Pajarito.jl.
• inactives::Vector = [] The indexes of inequality constraints that should be turned off.

MakeInputDict(m::Model)

Returns input_dict::Dict which contains the state space representation of the system equations stored in the flexibility model. The resulting input_dict is used by most of the flexibility analysis functions.

Arguments

• m::Model The flexibility model.

AddSystemExpressions(m::Model, input_dict::Dict, [num_scenarios::Int = 0])

Returns a vector of vectors where the first contains all the inequality expressions corresponding to the inequalities defined in input_dict and the second contains all of the equality expressions corresponding to the equalities defined in input_dict.

Arguments

• m::Model The flexibility model.
• input_dict::Dict Input dictionary as defined by MakeInputDict.
• num_scenarios::Int = 0 The number of scenerio subproblems, 0 turns off this feature.

RandomVariable(m::Model, mean::Number, name::AbstractString)

Return a RandomVariable DataType for given flexibility model m given the mean and the name. An anonymous JuMP variable is added directly to the flexibility model and its is appended to FlexibilityData.RVcols.

Arguments

• m::Model The flexibility model.
• mean::Number The variable mean.
• name::AbstractString The variable name.

RecourseVariable(m::Model, name::AbstractString)

Return a RecourseVariable DataType for given flexibility model m given the name. An anonymous JuMP variable is added directly to the flexibility model and its is appended to FlexibilityData.recourse_cols.

Arguments

• m::Model The flexibility model.
• name::AbstractString The variable name.

Base.show(io::IO, a::FlexibilityExpr)

Extend Base.show to print flexibility expressions.

JuMP.addconstraint(m::Model, constr::FlexibilityConstraint)

Extend the JuMP.addconstraint function to handle FlexibilityConstraint types.

JuMP.show(io::IO,c::FlexibilityConstraint)

Extend the JuMP.show function to handle FlexibilityConstraint types.

JuMP.constructconstraint!(flex_aff::FlexibilityExpr, sense::Symbol)

Extends JuMP.constructconstraint! for FlexibilityExpr types.

ComputeCenter(m::Model, center::Symbol, solver, toler::Number, only_positive::Bool)

Returns a center point that can be used to replace the mean if desired.

Arguments

• m::Model The flexibility model.
• center::Symbol Indicates the type of center, accepted arguments are :feasible and :analytic.
• solver The solver which must be an NLP solver for the analytic center.
• toler::Number The tolerance to check solution validity.
• only_positive::Bool Indicates if the center need by strictly positive.