Getting Started

This tutorial provides a step-by-step guide to modeling and solving a bilevel mixed-integer programming (BiMIP) problem using PowerBiMIP. By the end of this guide, you’ll understand the basic workflow and be able to apply it to your own optimization problems.

Installation

Before we begin, please make sure you have successfully installed PowerBiMIP and its dependencies (MATLAB, YALMIP, and a MILP solver like Gurobi). If you haven’t, please follow the instructions on the Installation page.

A Simple Bilevel Programming Example

Let’s start with a classic textbook example of a BiMIP problem. This will help us illustrate the core components of a PowerBiMIP model. The full script for this example can be found in examples/BiMIP_benchmarks/BiMIP_toy_example1.m.

1. Mathematical Formulation

The problem is composed of two levels: an upper-level problem and a lower-level problem. The decision variables of the upper level affect the constraints of the lower level, and the solution of the lower level, in turn, influences the upper level’s objective function.

Upper-Level Problem: The leader’s goal is to minimize its objective function by choosing the variable x.

min_{x}   -x - 10*z
s.t.
      x >= 0
      -25*x + 20*z <= 30
      x   + 2*z  <= 10
      2*x - z    <= 15
      2*x + 10*z >= 15

Lower-Level Problem: The follower’s problem is to minimize its objective by choosing variables y and z, given the upper-level decision x.

min_{y,z}   z + 1000 * sum(y)
s.t.
      -25*x + 20*z <= 30 + y(1)
      x   + 2*z  <= 10 + y(2)
      2*x - z    <= 15 + y(3)
      2*x + 10*z >= 15 - y(4)
      z >= 0
      y >= 0

2. MATLAB Implementation with PowerBiMIP

Now, let’s translate this mathematical model into MATLAB code using PowerBiMIP.

Step 1: Initialization

It’s good practice to start with a clean slate. This command clears the MATLAB workspace, closes all figures, clears the command window, and resets YALMIP’s internal memory.

clear; close all; clc; yalmip('clear');

Step 2: Define Variables

We use YALMIP’s syntax to define the decision variables. We organize them into a model structure with specific fields: var_upper for upper-level variables and var_lower for lower-level variables.

model.var_upper.x = intvar(1,1,'full'); % Upper-level integer variable
model.var_lower.z = intvar(1,1,'full'); % Lower-level integer variable
model.var_lower.y = sdpvar(4,1,'full'); % Lower-level continuous variables

Step 3: Formulate the Model

Next, we define the constraints and objective functions for both levels. The constraints are created as standard YALMIP objects. The model structure must strictly follow these field names: cons_upper, cons_lower, obj_upper, and obj_lower.

Upper-Level Constraints & Objective:

% --- Upper-Level Constraints ---
model.cons_upper = [];
model.cons_upper = model.cons_upper + (model.var_upper.x >= 0);
model.cons_upper = model.cons_upper + (-25 * model.var_upper.x + 20 * model.var_lower.z <= 30);
model.cons_upper = model.cons_upper + (model.var_upper.x + 2 * model.var_lower.z <= 10);
model.cons_upper = model.cons_upper + (2 * model.var_upper.x - model.var_lower.z <= 15);
model.cons_upper = model.cons_upper + (2 * model.var_upper.x + 10 * model.var_lower.z >= 15);

% --- Upper-Level Objective ---
model.obj_upper = -model.var_upper.x - 10 * model.var_lower.z;

Lower-Level Constraints & Objective:

% --- Lower-Level Constraints ---
model.cons_lower = [];
model.cons_lower = model.cons_lower + (-25 * model.var_upper.x + 20 * model.var_lower.z <= 30 + model.var_lower.y(1,1) );
model.cons_lower = model.cons_lower + (model.var_upper.x + 2 * model.var_lower.z <= 10 + model.var_lower.y(2,1) );
model.cons_lower = model.cons_lower + (2 * model.var_upper.x - model.var_lower.z <= 15 + model.var_lower.y(3,1) );
model.cons_lower = model.cons_lower + (2 * model.var_upper.x + 10 * model.var_lower.z >= 15 - model.var_lower.y(4,1) );
model.cons_lower = model.cons_lower + (model.var_lower.z >= 0);
model.cons_lower = model.cons_lower + (model.var_lower.y >= 0);

% --- Lower-Level Objective ---
model.obj_lower = model.var_lower.z + 1e3 * sum(model.var_lower.y,'all');

Step 4: Configure and Run the Solver

Finally, we configure the solver settings using BiMIPsettings and then call the main solver function solve_BiMIP. Note that the input to solve_BiMIP is now just the model structure and the ops options.

perspective: Set to 'optimistic' for this problem.
method: We use the 'exact_KKT' method, which is an exact algorithm.
solver: We specify 'gurobi' as our underlying MILP solver.
verbose: A setting of 2 provides detailed solver output.

% Configure solver settings
ops = BiMIPsettings( ...
    'perspective', 'optimistic', ...
    'method', 'exact_KKT', ...
    'solver', 'gurobi', ...
    'verbose', 2 ...
    );

% Call the solver
[Solution, BiMIP_record] = solve_BiMIP(model, ops);

3. Understanding the Results

After running the script, PowerBiMIP will return two main outputs:

Solution: A structure containing the optimal objective values.
BiMIP_record: A structure that records detailed information about the solution process, such as solving time and algorithm iterations.

For this example, the optimal solution is x = 2, z = 2, resulting in an upper-level objective value of -22.

Next Steps

Congratulations! 🎉 You’ve successfully solved your first bilevel problem with PowerBiMIP.

From here, you can:

Explore other examples in the examples/ folder.
Learn about Two-Stage Robust Optimization (TRO) features.
Apply PowerBiMIP to your own research problems.