Learning for Constrained Optimization: Identifying Optimal Active Constraint Sets
Sidhant Misra Line Roald Yeesian Ng













































1]Sidhant Misra
2]Line Roald
3]Yeesian Ng
[1]sidhant@lanl.gov
[2]roald@wisc.edu
[3]yeesian@mit.edu


















 
Learning for Constrained Optimization: 
 
Identifying Optimal Active Constraint Sets
    [
    December 30, 2023
====================================================================================




In many engineered systems, optimization is used for decision making at time-scales ranging from real-time operation to long-term planning. This process often involves solving similar optimization problems over and over again with slightly modified input parameters, often under tight latency requirements. 
We consider the problem of using the information available through this repeated solution process to directly learn a model of the optimal solution as a function of the input parameters, thus reducing the need to solve computationally expensive large-scale parametric programs in real time. Our proposed method is based on learning relevant sets of active constraints, from which the optimal solution can be obtained efficiently. Using active sets as features preserves information about the physics of the system, enables interpretable models, accounts for relevant safety constraints, and is easy to represent and encode.
However, the total number of active sets is also very large, as it grows exponentially with system size. 
The key contribution of this paper is a streaming algorithm that learns the relevant active sets from training samples consisting of the input parameters and the corresponding optimal solution, without any assumptions on the problem structure. 
The algorithm comes with theoretical performance guarantees, and is known to converge fast for problem instances with a small number of relevant active sets. It can thus be used to establish the practicability of the learning method.
Through extensive experiments on the Optimal Power Flow problem, we observe that often only a few active sets are relevant in practice, suggesting that the active sets is the appropriate level of abstraction for a learning algorithm to target.







































§ INTRODUCTION



The use of optimization methods to improve economic efficiency and enhance system performance is growing rapidly in engineering applications. Examples include the operation and planning of infrastructure like electric power grids, gas transmission systems, water and heating networks and transportation networks. Many of these systems operate in highly variable conditions, due to developments such as renewable energy integration in power systems and uncertain demand-side behaviour. To maintain safe and efficient operations, new optimized set-points are obtained by resolving the optimization problems as system conditions (and corresponding system parameters) change. The optimization problems exhibit two key features: 

    
  * There exists a mathematical model of the physical system, including objective function and technical constraints. For optimization purposes, this model is assumed to be fully known except for those input parameters whose values are obtained in real-time, giving rise to parametric programming problems. 
    
  * Constraint satisfaction corresponds to enforcing operational safety limits and adhering to physical laws, and is of paramount importance.

As the optimization problems are solved over and over again in response to changing conditions, we build a rich history of input parameters and their corresponding optimal solutions. 
However, solving these parametric programs in real-time for realistically sized systems can be prohibitively difficult, owing to the combination of complex physical laws and tight latency requirements. This limits the frequency with which the system set-points can be updated. 










The goal of this paper is to develop a framework that uses machine learning methods to enable these parametric programs to be solved more efficiently. 
An intuitive approach to apply machine learning would be to use data from existing solutions, and directly learn the optimal solution from the input data, as illustrated in Figure <ref>. Previous literature on representing the mapping from parameters to optimal solutions in parametric programs often used reinforcement learning (RL), where the underlying physical constraints of the problem are assumed to be unknown.
We refer the reader to <cit.> for examples of this approach using neural networks.
However, the neural networks mostly provide only an approximate model of the physical constraints, which may result in solutions that violate important system constraints. This has prompted a literature on “safety problems" <cit.>. Approaches to mitigate those safety concerns include modifying the exploration procedure by restricting it to a set of safe policies <cit.>, re-projecting the decisions made by the neural network into a known set of safe decisions <cit.>, or include an outer safety loop which overrides potentially unsafe decisions by the neural net <cit.>. Nonetheless, it remains difficult to enforce guarantees on behavioral constraints in the resulting policy <cit.> and verify that the resulting neural network policy provides a safe policy <cit.>. Approaches that restrict or limit the neural net also may produce sub-optimal solutions. 

A main drawback of these existing learning methods when applied in the context of engineered systems is hence their inability to enforce the constraints in the optimization accurately. One important reason for this drawback is that the methods are unable to leverage pre-existing knowledge about the mathematical form of the optimization problem. For best results, the learning methods used should exploit any known structure of the optimization problem in consideration. 












In this paper, we take a different approach to learning the optimal solution. Instead of directly learning the mapping from the input data to the optimal solution, we introduce the active set as an intermediate step in our learning procedure. To the best of our knowledge, this has never been explored before in the literature.
The optimal active set in an optimization problem corresponds to the set of constraints that are satisfied with equality in the corresponding optimal solution. Learning the optimal active set arguably constitute the right level of abstraction for applying machine learning algorithm to optimization problems for several reasons:

    
  (i) From an optimization perspective, the optimal active set in an optimization problem constitutes the minimal information required to recover the optimal solution. All optimization algorithms such as interior point methods, active set methods, and the simplex algorithm can see significant speed-up from the knowledge of the active set at optimality. In some cases, such as in linear programming (LP), the knowledge of the optimal active set can simplify the problem to the extent that it can be solved analytically, or in a single iteration.
    
  (ii)  From a machine learning perspective, using the active sets as features allows us to reduce the dimension of the learning task. Instead of learning a complex multi-dimensional, continuous-valued mapping from the input parameters to the optimal solution, the problem is converted to the simpler task of learning the mapping between the input parameters and a finite number of active sets, allowing us to exploit the known mathematical model of the system. Thus, in contrast to the standard approach in reinforcement learning, we spend no effort on learning the mathematical model itself.
    
  (iii) From an application perspective, active sets often correspond to meaningful operational patterns, allowing for interpretable models that can assist humans in decision-making. This point is critical in some engineering applications, where owing to the high stakes involved (such as avoiding black out in power grid), the final decision making is frequently performed by experienced system operators.
    
  (iv) Our intuition tells us that human engineered systems often admits just a few modes of operation. This is borne out in the numerical experiments that we ran on the optimal power flow (OPF) problem across a wide variety of transmission networks, which translates into only a small number of active sets being relevant in practice. This offers significant advantage for a machine learning approach to learn these active sets from data.





Therefore, a key contribution of this paper is to show that learning the active set is an efficient and viable approach in many practical problems, with low number of relevant active sets. To this end, we propose a statistical learning approach which provides probabilistic guarantees for out-of-sample performance, and characterizes the complexity of the task in terms of the number of active sets required to produce an optimal solution with high probability. 
The input to the algorithm is a set of training samples consisting of the input parameters and the corresponding optimal solution, without any assumptions on the problem structure. 
The algorithm is designed to terminate when a sufficient number of optimal active sets have been discovered. The algorithm also allows for the user to set termination criteria that provides theoretical performance guarantees. In the worst case, the algorithm may require an exponential number of samples. However, we prove that it will terminate fast for so-called low-complexity systems where the number of optimal active sets is small. 
The algorithm thus simultaneously identifies the optimal active sets and establishes the practicability of the active-set-based learning method, which hinges on only a modest number of optimal active sets.
As soon as the relevant active sets are identified, we can use them to predict the optimal solution, as outlined in Figure <ref>.








A related, but distinct field of study is explicit Model Predictive Control (MPC), which was developed to improve solution times of MPC problems in an online setting with time-varying parameters. For parametric programs with a convex quadratic objective and linear constraints, the parameter space can be partitioned into a number of regions (corresponding to different optimal active sets), where the optimal solution in each region is an affine function of the input parameters.
Explicit MPC approaches generally involve an offline procedure to identify the full set of regions and corresponding optimal solutions <cit.>. They can then be stored in sophisticated data-structures to reduce memory storage requirements <cit.>. There are also advances in fast search algorithms for the “point location problem" of identifying the corresponding active set for a given state point based on binary search trees <cit.> and other data structures <cit.>. 
Further, <cit.> propose to design a deep neural net to specifically captures the piece-wise affine structure of the MPC solution manifold, and provides theoretical bounds on the number of neurons and layers to enable an exact representation of the solution.

An important difference to our approach is that explicit MPC searches through the full parameter space, and is hence limited to low dimensional systems, whereas our method uses information about the distribution over the input parameters to only investigate the practically relevant parts of the input space. 
Further, our approach is more general – explicit MPC is focused on solving convex quadratic programs, whereas the type of problems under consideration in this paper belong to a more general class of problems, sometimes including non-convex or integer constraints. 



In the following we lay down the foundations for a statistical learning approach which provides probabilistic guarantees for out-of-sample performance, and characterizes the complexity of the task in terms of the number of active sets required to produce an optimal solution with high probability. 
Our main contributions in this paper are as follows:

(i) Establishing active sets as a natural and effective means to learn solutions to general parametric optimization problems. To the best of our knowledge, this has never been explored before in the literature.

(ii) Developing a streaming algorithm termed DiscoverMass with rigorous performance guarantees to discover the collection of relevant active sets, without any a priori assumptions on problem structure. 
The algorithm serves a dual purpose, as it both learns the relevant active sets and provides evidence of whether the suggested active-set-based learning method is applicable to a given problem instance.

(iii) Establishing the viability of the active-set approach on a practical problem. We employ DiscoverMass on a set of Optimal Power Flow benchmarks representative of optimization problems solved in operation of electric transmission grids to demonstrate that the number of relevant active sets is indeed small and learning them is highly efficient in terms of the number of samples required.

We emphasize that this paper does not provide a full end-to-end learning framework based on learning the active set. While we provide examples for how the discovered collection of active sets can be used to recover the optimal solution for a new realization of the input parameters using an ensemble policy or classification, the development of those frameworks will be considered in future work. In this paper, we lay the foundation for this future work by showing that methods based on learning the active sets are viable for practical problems. 

The remainder of the paper is organized as follows. We first describe the main idea behind learning active sets in Section <ref>. We then develop a learning algorithm to identify a collection of active sets from training samples in Section <ref> and establish theoretical performance guarantees. We demonstrate the effectiveness of our method using an application to power systems control in  Section <ref>. Section <ref> concludes the paper.



















































































































































§ INTRODUCTION





In this paper, we are interested in learning interpret-able representations for control and planning. We begin from the standpoint of being able to solve optimization problems corresponding to convex parametric programs based on historical realizations of the input parameters, and pose the question: what can we learn from the process?

This question arises as the use of optimization methods to improve economic efficiency and enhance system performance is growing rapidly in many practical engineered systems, particularly in infrastructure networks such as the electric power grid, gas transmission systems, water and heating networks and transportation networks. Many of these systems must be operated in the face of increasing uncertainty, due to developments such as renewable energy integration in power systems, imperfect system information due to aging components, and uncertain demand-side behaviour. However, solving these operational planning problems in real-time systems can be prohibitively difficult, as they often need to be solved for large systems in a very limited time.








 §.§ Literature Review

With sustained developments in optimization technology and increasing availability of powerful optimization solvers such as Ipopt <cit.>, KNITRO <cit.>, Gurobi <cit.>, MOSEK <cit.>, SCIP <cit.>, OSQP <cit.> or ECOS <cit.>, the latency of solving large-scale optimization problems has decreased tremendously <cit.>. 


This makes them practical for solving offline strategic planning problems. However, these solvers often use active-set or interior-point methods that are sequential in nature, making them difficult to effectively parallelize for real-time operational decision planning.

This prompts two broad research directions with the goal of representing the mapping from parameters to optimal solutions in convex parametric programs: the first is to assume that the underlying physical constraints of the problem are unknown, and to learn it through interactions with an environment in reinforcement learning (RL). The second approach is to assume that the constraints of the problem are fully specified by a parametric program in Model Predictive Control (MPC), and elicit or approximate it using known structure in the underlying geometry.


Recent advances in RL involves training neuralnets to either approximate the value function or the policy network <cit.>. However, the neuralnet might not reflect the physical constraints accurately, and can result in infeasible solutions due to violations of safety constraints, prompting a literature on safety problems that can arise in AI planning <cit.>. Some approaches to mitigate it include modifying the exploration procedure to restrict it to a set of safe policies <cit.>, or re-projecting from potentially unsafe decisions made by the neuralnet policy into a known set of safe decisions <cit.>. Nonetheless, it remains difficult to enforce guarantees on behavioral constraints in the resulting policy <cit.> and verify that the resulting neuralnet policy corresponds to a safe policy <cit.>. More fundamentally, it does not sidestep the need for algorithms in RL to explore and interact with an environment to observe the consequences of different decisions and alter their behavior accordingly <cit.>.






In Explicit MPC, techniques have been developed to elicit the mapping from parameters to optimal solutions in parametric programs <cit.>, allowing them to be programmed into embedded devices. They typically rely on efficient methods for enumerating over all possible bases, with implementations in software such as the Multi-Parametric Toolbox <cit.>. However, due to the combinatorial explosion in possible configurations for the bases, they are often only computationally feasible for small systems. There are ongoing efforts to learn an approximate policy <cit.>, with applications in autonomous driving <cit.> among others, but there has been little connection established between the techniques developed there and the notions of basis functions and generalization in machine learning.



 §.§ Our Contributions


In this paper, we study a data-driven approach to the learning the mapping between uncertain parameters and the optimal solution in convex parametric programs, where we incorporate the safety constraints directly into the model.
The key observation in this approach is to define the basis functions based on the active set of constraints, i.e., each active set corresponds to a basis function.
In addition to easily accounting for the safety constraints, the active sets, by nature of their construction often correspond to meaningful systems states for domain experts to interpret <cit.>. 
In essence, we effectively split the challenging learning task into simpler and interpret-able steps consisting of first learning the set of active sets relevant to practice, then learning the mapping from the uncertain parameter to the corresponding active set and finally the mapping from the active set to the optimal solution. 







The rest of this paper is organized as follows: we describe how to recover an optimal solution from a collection of active sets in section <ref>, before adopting a statistical learning approach to elicit a sparse set of active sets that would approximate the underlying decision-model in section <ref>. We provide probabilistic guarantees on its generalization capabilities in section <ref>, before demonstrating its application to power systems control in section <ref>.














Active Set Learning for Optimization 

The aim of this paper is to learn decision policies x^*(ω) which returns the optimal solution x^* for parametric programs of the form

    x^*(ω)∈x∈(ω)argmin f(x,ω)

where x denotes the decision variables, ω represent the uncertain input parameters, and

    (ω)={x∈^n: g_j(x,ω)≤ 0   ∀ j∈ [m]}

defines the set of feasible solutions, and f,g_1,…,g_m are functions in x for all ω. We restrict ourselves to a set of parameters ω∈Ω, for which the optimization problem is actually feasible:

    Ω={ω:𝕏(ω)≠∅}.

For any given solution x∈𝕏(ω), we define the active set corresponding to x as the set of constraints  that are binding at x.
For a given parameter realization ω, we define the optimal active set ^*(ω) as the active set at the optimal solution x^*(ω) [When there are multiple optimal solutions, the degeneracy can be handled by using any consistent tie-breaking rule,
such that the same realization ω will always be mapped to the same optimal active set.]


If we know the optimal set ^* for a given ω∈Ω, we can solve the following optimization problem

    x^*_^*(ω)∈x∈𝕏_^*(ω)argmin  f(x,ω)

based on the reduced set of constraints ^*

    𝕏_^*(ω) = {x∈^n: g_j(x,ω) ≤ 0   ∀ j∈^*}

The solution to (<ref>) will be feasible and optimal for the original problem (<ref>), but is easier to obtain since we typically have that |^*| ≪ m. 
If we have a method to map the uncertain parameters ω to the corresponding optimal active set ^*, we are thus able to obtain a the optimal solution efficiently.

Consider the case where we draw a new sample ω and do not know the optimal active set ^*, but have a collection of possible candidates _1,...,_N. Then, it is possible to obtain the optimal solution by solving the reduced problem (<ref>) for each active set and check the resulting solutions for feasibility. Since each reduced problem (<ref>) is a relaxation of (<ref>), the solutions will either be optimal or infeasible for the original problem. 


































 §.§ Learning the collection of relevant active sets

While the collection of all possible active sets is finite, it is also exponential in the problem size, making it a potentially complex learning task with high sample complexity. However, in many applications such as infrastructure systems that were engineered to accommodate particular modes of operation, only a few active sets are relevant in practice. The learning task at hand thus becomes to identify those important active sets.

In our approach, the important active sets are the active sets that have the highest probability of being optimal, under a given probability distribution over the uncertain input parameters ω. Our procedure to discover the active sets is illustrated in Figure <ref>. We first draw i.i.d. training samples ω_i from the probability distribution over ω. For each sample, we solve the optimization problem (<ref>) to obtain the corresponding optimal solution and active set _i, which is now marked as an observed active set (in color in Figure <ref>). The collection of all active sets _1, ..., _N that were observed during the first ω_i, ...,ω_N samples form the set of observed active sets. The still unobserved active sets, marked in grey in Figure <ref>, are the active sets which were not optimal for any of our i.i.d. samples ω_i.



 §.§ From active sets to the optimal solution

The knowledge of the relevant active sets can be used to more efficiently obtain the optimal solution. We describe two possible approaches.



  
Ensemble policy
The ensemble policy obtains the optimal solution by solving the reduced problem (<ref>) in parallel for each relevant active set , and then verifying the resulting solution for feasibility in (<ref>). This method is described in detail in <cit.>, and illustrated in Figure <ref>. 



  
Classification
This approach is based on training a classifier to learn the mapping from the parameter ω to the corresponding optimal active set. The number of classes in this classifier is same as the number of relevant active sets. Once the optimal active set is obtained, the optimal solution can be obtained by solving the reduced optimization problem in (<ref>). 

Both the above approaches are simpler and more efficient when the number of potential active sets is small.
In this paper we utilize the ensemble policy in our experiments to test the performance of our method, and leave classification as a topic for future research.

We envision any learning approach based on learning active sets to rely on combining a new realization  ω with one or more observed active sets, as illustrated in Figure <ref>. 
If the optimal active set corresponds to one of the observed active sets, then the optimal solution lies within the set of optimal solutions to the corresponding collection of relaxed problems (<ref>). On the other hand, if the optimal active set of the new sample was previously unobserved, the solutions to the all relaxed problems (<ref>) with _i,...,_N may be infeasible. Since this is an undesirable outcome, we would like to bound the probability of these outcomes. 
A main contribution of this paper is to provide a learning algorithm which guarantees that the probability of encountering a previously unobserved optimal active set falls below a pre-specified level. We develop a streaming algorithm which draws new samples and observes new optimal active sets until the desired performance level can be guaranteed. This helps us distinguish between cases where we have discovered active sets which cover a smaller or larger portion of the feasible space, as illustrated in Figure <ref>. The algorithm is formalized in Section <ref>.





















































































































Algorithm to learn relevant active sets 
In this section, we describe a streaming algorithm that enables us to learn a collection of active sets whose combined probability mass exceeds a user-defined safety level. 
The algorithm draws random samples of the uncertain input parameters ω_i from the distribution _ω(.), solves the optimization problem (<ref>) to obtain the corresponding active set _i and uses a stopping criterion to ensure that sufficiently many active sets have been discovered. We provide details of the algorithm, and theoretical guarantees on its performance below.


Let ℬ denote the finite set of all possible active sets, which is exponential in the number of constraints in the problem. 
We define observed and unobserved active sets in the following way. 
 
Let ω_1, …, ω_M be M i.i.d. samples drawn from the uncertainty distribution and _1, …, _M denote the optimal active sets corresponding to each ω_i. We call 𝒪_M = ∪_i=1^M {_i} the set of observed active sets, and 𝒰_M = ℬ∖𝒪_M the set of unobserved active sets.

In Figure <ref>, 𝒪_M and 𝒰_M correspond to the colored and white regions respectively. The streaming algorithm iteratively increases the number of samples M until the observed set 𝒪_M contains most of the important active sets. A crucial aspect of the algorithm is that knowledge of neither the shape and number of the active sets (colored regions) nor the number of samples M is required a-priori, but is decided on-the-fly. The algorithm therefore is guaranteed to converge for any system if enough samples are available, but can also be stopped at any time before termination with a lower safety guarantee.

In the following, we use the notion that the importance of an active set or a collection of active sets corresponds to the mass contained in it.
 
The mass of an active set is defined as the probability that it is optimal for (<ref>) under the distribution _ω(.) on ω. For ∈ℬ,

    π() = _ω( = ^*(ω))

For any subset 𝒮⊆ℬ, we define the mass of 𝒮 as π(𝒮) = ∑_∈𝒮π().

After drawing M samples and observing a collection of active sets 𝒪_M, we compute the rate of discovery. The rate of discovery quantifies the fraction of samples that correspond to observing an active set that was not observed among the first M samples, and is formally defined below.
 
Let W be a positive integer denoting the window size. Let ω_1, …, ω_M+W be M+W i.i.d. samples drawn from the uncertainty distribution and let _i denote the optimal active set for ω_i. We denote by X_i the random variable that encodes whether a “new active set", i.e., an active set not observed within the first M samples, was observed in the (M+i)^th sample, i.e.
 
    X_i = 
                        1,   _M+i∉{_1}∪…∪{_M}, 
    
                        0,  .

    Then the rate of discovery over the window of size W is given by ℛ_M,W and is defined as
    
    ℛ_M,W  = 1/W∑_i=1^W X_i.



The rate of discovery ℛ_M,W is directly related to the mass of the undiscovered set 𝒰_M of active sets. Intuitively, if we have already discovered the active sets that contain most of the mass, then the rate of discovery is expected to be low since most of the time we are likely to observe one of the active sets we have already observed so far. The following theorem makes this notion rigorous.
 
Let {ω_i}_i=1^∞ be a sequence of i.i.d. observations and define the unobserved set and rate of discovery as in Definition <ref> and <ref> respectively. Then the following statements hold.
Let the probability mass in the unobserved set be denoted by π(𝒰_M) and let the window size be such that 
    
    W = W_M ≥ c max{logM,  log M}.

    Then we can bound the probability mass in the undiscovered set by
    
    ℙ   (π(𝒰_M) - ℛ_M,W≤ϵ  ∀ M ≥ 1 ) > 1 - γ/γ-11/(M-1)^γ-1,   γ = c ϵ^2/2.


The proof is deferred to Section <ref>.
Theorem <ref> forms the basis of the stopping criterion of our streaming algorithm. In the following subsection, we present this algorithm and provide rigorous guarantees on its performance based on the rate of discovery criterion. 


Streaming Algorithm for Learning Important Active Sets 
We present the following algorithm called DiscoverMass which takes in as input parameters the quantities α, ϵ, δ >0. Here, α represents the limit on the undiscovered mass and requires that the algorithm returns a set of active sets of mass at least α. The quantity ϵ represents the difference between the undiscovered mass and the rate of discovery, and 1-δ represents the confidence.




We can also choose the hyperparameter γ. We have explicitly left the dependence on all the parameters in the algorithm description so that interested readers can choose the values that may be suitable for their application. We specify the concrete set of chosen values along with a brief analysis of the effect of tuning the hyperparameters in the numerical experiment section.
 
    
DiscoverMass terminates with a set of observed active sets 𝒪 that has at least 1-α of the mass with high probability. The following theorem guarantees that with the stopping criterion described above, DiscoverMass indeed succeeds in meeting the requirements imposed by the input parameters. 

 
If the algorithm DiscoverMass is employed with input parameters α, ϵ, δ, γ, then with probability at least 1-δ the algorithm will terminate having discovered active sets with mass at least π(𝒪) ≥ 1 - α.

The proof of Theorem <ref> can be found in Section <ref>.
A key advantage of DiscoverMass supported by Theorem <ref> is that it makes absolutely no a priori assumption regarding the underlying distribution of active sets in the system. This makes DiscoverMass a rigorous tool to quantify the unknown underlying distribution.


    
    
    
    
    


Low-complexity systems 
Although Theorem <ref> guarantees that the algorithm DiscoverMass will almost always succeed in finding the desired amount of mass, it gives no guarantees on the number of steps M the algorithm needs to discover this mass. In fact, the number of steps is highly dependent on the properties of the system and can vary significantly across systems. For example in a system where there are 2^N active sets all of which are equally probable, i.e. _ω(_i) = 2^-N for all i ∈ [2^N], then it is clear that discovering (1-α) fraction of the mass requires at least M > (1-α)2^N samples. However, for some systems, which we refer to as low-complexity systems, it is possible to discover a large fraction of the mass quickly. Fortunately, many practical systems fall within this category. In Section <ref>, we have found that many practical systems can be described as low-complexity as it is possible to discover a large fraction of the mass quickly. We quantify this intuition with the following definition.
    
     
        We say that a system is a low-complexity system with parameters α_0 and K_0 if there are at most K_0 active sets that contain at least 1-α_0 fraction of the mass. More precisely, there exists active sets _i_1, …, _i_K_0 such that ∑_j=1^K_0π(_i_j) > 1-α_0.
    
    
    The next theorem guarantees that for low-complexity systems the algorithm terminates within a few iterations.
     
        Suppose that we run the algorithm DiscoverMass(α,ϵ,δ,γ) on a low-complexity system with parameters (α_0,K_0) with α >  α_0. Then with probability at least 1 - δ - δ_0 the algorithm terminates in less iterations than M = 1/α - α_0(K_0 log 2 + log 1/δ_0 ).
    
    The proof of Theorem <ref> is found in the next section. From Theorem <ref>, we see that the sample complexity of discovering active sets with sufficient mass for low-complexity systems scales at most linearly with the number of underlying important active sets K_0. 
    
  


Proofs 
In this section, we provide the proofs of the three theorems stated in Section <ref>.


 §.§ Proof of Theorem <ref>

We first prove part (a).  For any M we can bound the probability that the rate of discovery is far from the mass of the unobserved set by 
 
    ℙ    (π(𝒰_M) - ℛ_M,W  >  ϵ ) 
       = ∑_uℙ  (π(𝒰_M) - ℛ_M,W  >  ϵ|𝒰_M = u)  ℙ(𝒰_M = u) 
       = ∑_uℙ(  π(u) -  ∑_i=1^W_M X_i /W_M>  ϵ|𝒰_M = u)  ℙ(𝒰_M = u),

        where the random variables X_i are defined in (<ref>) and conditioned on 𝒰_M = u are i.i.d. Bernoulli random variables with mean π(u). Using standard large deviation inequalities, we can bound the probability in (<ref>) as
        
    ℙ  (π(𝒰_M) - ℛ_M,W  >  ϵ )     < ∑_u e^-W_Mϵ^2/2ℙ(𝒰_M = u)  
       = e^-W_Mϵ^2/2.

        Therefore,
        
    ℙ    (∃ M  ≥ 1      π(𝒰_M) - ℛ_M,W  >  ϵ ) 
       ≤∑_M=1^∞ℙ (π(𝒰_M) - ℛ_M,W  >  ϵ ) 
       = ∑_M=1^Mℙ (π(𝒰_M) - ℛ_M,W  >  ϵ ) 
       +  ∑_M=M^∞ℙ (π(𝒰_M) - ℛ_M,W  >  ϵ )  
        < ∑_M=1^M  e^-c logMϵ^2/2 + ∑_M=M^∞ e^-c log M ϵ^2/2
       ≤1/M^c ϵ^2/2 -1 + ∑_M=M^∞1/M^c ϵ^2/2≤γ/γ-11/(M-1)^γ-1.

       
        
        
        
           
        
   


 §.§ Proof of Theorem <ref>

    Let 𝔼 be the event that the algorithm terminates after having discovered smaller than ϵ fraction of the mass. Then
        
    ℙ( 𝔼)    ≤ℙ(∃ M ≥ 1 |ℛ_M,W < ϵ - δ, π(𝒰_M)  > ϵ) 
       ≤ℙ(∃ M ≥ 1 |π(𝒰_M) - ℛ_M,W  > ϵ) 
       (a)<γ/γ-11/(M-1)^γ-1(b)=δ,

        where the last inequality (a) follows from Theorem <ref>, and (b) follows from the value of M in the initialization step in DiscoverMass.
        


 §.§ Proof of Theorem <ref>

	    Let I^* = {_i_1, …, _i_K_0} denote the set of K_0 active sets defined in Definition <ref> that contain at least 1-α_0 of the mass.  For any M ≥ 1 recall that 𝒰_M denotes the set of unobserved active sets as in Definition <ref>. We define the sets  𝒰_M_I^* = 𝒰_M ∩ I^* and 𝒰_M_I̅^* = 𝒰_M ∩I̅^* where I̅^* denotes the complement of the ensemble I^*, i.e., the set of active sets that are not a part of I^*. It follows that
	    𝒰_M = 𝒰_M_I^*∪𝒰_M_I̅^* and π(𝒰_M) = π(𝒰_M_I^*)  + π(𝒰_M_I̅^*).
	    We compute
	    
    ℙ   (π(𝒪_ℳ) < 1 - α) = ℙ(π(𝒰_M) > α) 
        = ℙ(π(𝒰_M_I^*)  + π(𝒰_M_I̅^*) > α) 
       ≤ℙ(π(𝒰_M_I^*) > α-α_0) 
       = ∑_u ⊆ I^*ℙ(𝒰_M_I^* = u)    1_{π(u) > α-α_0}
       ≤∑_u ∈ I^*(1- (α-α_0) )^M 
       ≤ 2^K e^-(α-α_0)M(a)<δ_0,

	    where (a) follows from the choice of M in the theorem. The proof follows by combining (<ref>) with Theorem <ref>.




Numerical results 

We test our learning framework on the DC Optimal Power Flow (OPF) problem, an optimization problem widely used in electricity market clearing as well as operation and planning of electric transmission grids. 
For real time operation, the reaction of generators to fluctuations in renewable energy must be calculated reliably within a short period of time, so it is an example with high practical relevance that is subject to stringent limits on computational time. 
A simple version of the DC OPF can be stated as follows



 
	
    ρ(ω)∈p∈^nargmin    c^⊤ p 
    s.t.    e^⊤(p-d+ω) = 0 
       ≤ p≤
       ≤ M(Hp-d+ω)≤


Here, decisions p∈^n are made on the active power generation at each generator. The aim is to provide power at minimum generation cost c^⊤ p to a set of loads d∈^v with random real-time variations ω. 

The solution is subject to a total power balance constraint (<ref>), as well as limits both on the minimum and maximum power generation ,∈^n at each generator (<ref>), and the minimum and maximum admissible flow ,∈^m on transmission lines (<ref>). The topology of the network is encoded in the matrix H∈^v× n mapping the power from each generator to their corresponding bus, and the matrix of power transfer distribution factors M∈^m× v <cit.>. Here e is the vector of ones.



Test set-up
We test our algorithm by running extensive simulations across a variety of networks <cit.> from the IEEE PES PGLib-OPF v17.08 benchmark library <cit.>. We report general results for 15 different test cases, varying in size from 3 to 1951 buses.


For each system, we consider two different distributions for the uncertain load deviations ω.
First, we assume a multivariate normal distribution, where ω is defined as random vector of independent, zero mean variables with standard deviations σ=0.03d and zero correlation between loads. 

Second, we assume that each entry ω_i follows a uniform distribution with support ω_i∈[-3σ, +3σ]=[-0.09d,0.09d]. The uniform distribution is intentionally designed to spread the probability mass more evenly across a larger region, thus enabling us to compare cases with a different number of relevant active sets and different probability mass in each active set. Note that our method is able to handle any distribution, as long as a sufficient number of training samples is available.

The algorithm DiscoverMass has four input parameters, α, δ, ϵ, γ. We set the maximum mass of the undiscovered bases α=0.05 and the confidence level δ=0.01. 
The remaining parameters are set to ϵ=0.04 and γ=2.
Note that ϵ and γ are hyperparameters that can be tuned, leading to different minimum window size and M. For details about the choice of parameters and their effects, we refer the reader to supplementary material. 
We run the algorithm until termination, or until M=22'000. 


The procedure is implemented in Julia v0.6 <cit.>, using JuMP v0.17 <cit.>, PowerModels.jl v0.5 <cit.> and OPFRecourse.jl <cit.>.





	
    

	
    


 §.§ Parameter choices in experiments

Varying the hyper-parameter γ results in different values for the window size W_M and minimum M, while changing ϵ influences both the stopping criterion R_M,W and the size of the window W. We note that there is a trade-off in the choice of both of those parameters. If γ is chosen such that M is low, the initial window size is smaller, but will rapidly start increasing as soon as M>M. For ϵ a lower value implies a steep increase in the window size, but also reduces the requirements on the rate of discovery.

The difference between the probability mass of the undiscovered bases and the rate of discovery R_M,W in the window W are set to ϵ = 0.04, while we choose the hyperparameter γ=2. Together, this implies the stopping criterion R_M,W≤0.01. The minimum window size W_M≥ 13'259 samples, and remains constant as long as M≤M=201 samples. We run the algorithm until termination, or until M=22'000. To assess how accurately we capture the probability of the undiscovered active sets, we run an out-of-sample test with 20'000 samples. 

Numerical results

Table <ref> lists the simulation results for the 15 different networks with sizes ranging from 3 to 1951 nodes. For each network, we report the number of samples to termination M, the number of active sets K_M discovered within the first M samples and the window size W_M required to guarantee a confidence of α<0.05. 
In addition, we list the rate of discovery R_M,W at termination. All results are included for both the normal and uniform distributed parameters, and the system is sorted based on the number of active sets. 
Note that the choice of input parameters leads to a minimum window size W_M≥ 13'259 samples, and remains constant as long as M≤M=201 samples. 



  
Low-complexity systems We observe that for most systems, the learning algorithm performs well. Most systems have a relatively low number of relevant active sets (<10), and terminate after less than 200 samples while the window size is still W_M = W_M,min = 13'259. System size is not a determining factor for the number of relevant active sets, with even the largest systems case1888_rte and case1951_rte exhibiting a low number of relevant active sets.
This validates the intuition that many engineered systems (at least arising in the power systems benchmarks we've looked at) can be considered low-complexity as defined in Section <ref>. Comparing M and K_M in Table <ref>, we see a strong correlation between the number of active sets with majority of the mass and the samples till termination M. This is in agreement with the sample complexity derived in Theorem <ref>. The fact that most systems are low-complexity suggests that our proposed method for parameterizing and learning an optimal policy based on active sets might be effective in a wide variety of settings arising in engineered systems.



We observe that the primary factor affecting the number of samples till termination is not the number of important active sets, but the existence of active sets with small but not insignificant probability mass (typically in the range 0.001-0.01). Due to their low probability mass, those active sets require a large number of samples for discovery, while the cumulative probability of the still undiscovered active sets is still high enough to exceed the stopping criterion R_M,W≤0.01.







  
High-complexity systems The outliers in this analysis are cases such as case24_ieee_rts, case73_ieee_rts, case200_pserc and case240_pserc.
For case73_ieee_rts (uniform) and case240_pserc (normal and uniform), the algorithm does not terminate until we reach the upper limit of M=22'000. For case240_pserc, the rate of discovery remains above R_M,W≥0.07 even after the discovery of 2993 active sets. 
We believe that this performance deficiency is due to atypical system characteristics. The two systems case200_pserc and case240_pserc, are the result of power grid models that have been subjected to network reduction, hence altering the properties of the systems relative to the other, more typical systems. Similarly, case24_ieee_rts and case73_ieee_rts are special in that they have a large number of generators relative to the number of nodes, and case73_ieee_rts consists of three identical areas. Due to the nature of their construction, these instances are useful adversarial systems for analyzing the performance of different learning algorithms.











  
Influence of parameter distribution We observe that the uniform distribution typically discovers a larger number of active sets compared to the normal distribution and requires more samples M before termination. 
With the probability mass spread more evenly across the probability space, there are a larger number of active sets with significant probability mass, and fewer active sets with high probability mass (for the multivariate normal, there is a concentration of mass in active sets closer to the forecasted value). However, for systems like case24_ieee_rts, the uniform distribution appears to encourage earlier termination due to a faster exploration rate. 






























  
Out of sample prediction performance
The learning algorithm is compared against a full optimization approach where for each realization of ω the problem (<ref>) is solved using standard optimization methods. This approach has higher computational requirements, but produces optimal solutions with probability one. Our learning based method has better computational efficiency but has a non-zero probability of failure, corresponding to returning infeasible solutions.



To assess the performance of our learning algorithm we evaluated the predictions generated by the ensemble policy (where the solution is evaluated for all discovered active sets, then checked for feasibility) on 20'000 test samples. 

The results are shown in Table <ref> which shows the probability of failure. We see that this value is consistently lower than the mass of the undiscovered set α provided as input to the DiscoverMass algorithm.


Intuitively, this is to be expected, since Theorem <ref> strives to provide  upper bounds on the mass of the undiscovered set that are universal in nature and are guaranteed to work for all almost any continuous optimization problem and various kinds of uncertainty, without pre-specifying the number of samples M. 






	






Conclusion

We present an approach to learn the mapping from a set of uncertain input parameters to the optimal solution for parametric programs. Our target application is engineered systems, where optimization problems used in operational decision making are typically solved repeatedly, but with varying input parameters. Our approach is based on learning the active sets of constraints at the optimal solution. By relying on active sets, we are able to decompose the learning task into a simpler classification task, and at the same time exploiting the knowledge of the exact mathematical model of the optimization problem, and enabling tight constraint enforcement. The viability of learning active sets and the performance of our approach is validated using a series of experiments based on benchmark optimization problem used in electric grid operation.


The prediction step in this paper used the so-called ensemble policy. In future work we will explore classification and ranking algorithms to learn the mapping between the realization of the uncertain parameter onto the optimal active set. Further topics of exploration include application to long-term planning and multi-stage optimization problems. 




























    

    

    




































    

    

    
















    














    

    

    






    

    







    













	

	

	

	

































	

	






informs2014
















































