Monte Carlo Optimization for Solving
Multilevel Stackelberg Games

inst1]Pravesh Koirala

[inst1]organization=Vanderbilt University,

inst1]Forrest Laine

Stackelberg games originate where there are market leaders and followers, and the actions of leaders influence the behavior of the followers. Mathematical modelling of such games results in what's called a Bilevel Optimization problem. There is an entire area of research dedicated to analyzing and solving Bilevel Optimization problems which are often complex, and finding solutions for such problems is known to be NP-Hard. A generalization of Stackelberg games is a Multilevel Stackelberg game where we may have nested leaders and followers, such that a follower is, in turn, a leader for all lower-level players. These problems are much more difficult to solve, and existing solution approaches typically require extensive cooperation between the players (which generally can't be assumed) or make restrictive assumptions about the structure of the problem. In this paper, we present a stochastic algorithm to approximate the local equilibrium solutions for these Multilevel games. We then construct a few examples of such Multilevel problems, including: a) a nested toll-setting problem; and b) an adversarial initial condition determination problem for Robust Trajectory Optimization. We test our algorithm on our constructed problems as well as some trilevel problems from the literature, and show that it is able to approximate the optimum solutions for these problems within a reasonable error margin. We also provide an asymptotic proof for the convergence of the algorithm and empirically analyze its accuracy and convergence speed for different parameters. Lastly, we compare it with existing solution strategies from the literature and demonstrate that it outperforms them.

Stackelberg games Multilevel Optimization Monte-Carlo algorithm Trajectory optimization Adversarial optimization

§ INTRODUCTION

Stackelberg Equilibriums are well-known and extensively studied economic phenomena. In their most rudimentary form, they occur when there is a market leader whose decision influences one or many market followers. These leaders and followers are constrained in their own way and are assumed to be rational players who seek to minimize their costs (or maximize their profits) while satisfying their constraints. Mathematical modeling of these games gives rise to a Bilevel Optimization problem of the following form:

\begin{align*}
\min_{x_1\in \real^{n_1}, x_2 \in \real^{n_2}} ~~&f^1(x_1, x_2) \\
s.t. ~~&g^1(x_1, x_2) \ge 0\\
&x_2 \in \arg\min_{x_2 \in \real^{n_2}} ~~f^2(x_1, x_2) \\
&~~~~~~~~~~~~~~~~~s.t.~~g^2(x_1, x_2)\ge 0
\end{align*}

Where the upper level player with the objective $f^1(x_1, x_2): \real^{n_1+n_2} \mapsto \real$ and constraints $g^1(x_1, x_2): \real^{n_1+n_2}\mapsto \real^{m_1}$ optimizes over $x_1 \in \real^{n_1}$ knowing that its choice of $x_1$ causes the lower-level player to adapt its response variable $x_2\in \real^{n_2}$ to minimize its objective $f^2(x_1, x_2): \real^{n_1+n_2}\mapsto \real$ subject to its constraints $g^2(x_1, x_2):\real^{n_1+n_2}\mapsto\real^{m_2}$. It is also generally assumed that the objective functions $f^1$ and $f^2$ and the constraints $g^1$ and $g^2$ are twice differentiable. But even with these assumptions, solution set for problems of this form generate not only non-convex but also non-smooth manifolds. In fact, finding an equilibrium point or a solution to these problems is known to be NP-Hard [Ben-Ayed and Blair, 1990, Blair, 1992]. Popular strategies to solve these problems are Vertex Enumeration methods [Bialas and Karwan, 1984], Complementary Pivoting methods [Júdice and Faustino, 1992], Mixed Integer Programming or Branch and Bound methods [Bard and Moore, 1990], and meta-heuristics based methods such as Genetic Algorithms [Oduguwa and Roy, 2002] and Particle Swarm optimization [Han et al., 2016] etc.

The problem discussed above is called a Bilevel problem because it has two levels of optimizers (alternatively referred to as decision makers or players in this text) with their own sets of decision variables and constraints. A natural extension of such a leader-follower game is, then, a Multilevel Stackelberg game that can be modeled as:

\begin{align*}
&Level_1& &~~~~~~~~~~~\hdots\\
&& &~~~~~~~~~~~~~\vdots\\
&Level_l & \min_{x_l, x_{l+1}..., x_L} ~~&f^l(X) \\
&& s.t. ~~&g^l(X) \ge 0\\
&Level_{l+1} & &x_{l+1},...,x_L \in \arg\min_{x_{l+1},...,x_L} ~~f^{(l+1)}(X) \\
&& &~~~~~~~~~~~~~~~~~~~~s.t.~~g^{(l+1)}(X)\ge 0\\
&& &~~~~~~~~~~~~~\vdots\\
&Level_{L}& &~~~~~~~~~~~\hdots\\
\end{align*}

Where $l \in {1, 2, ..., L}$ indicates player level and $x_l \in \real^{n_l}$ is the variable that player $l$ controls. Similarly, $X \in \real^{n}$ (where $n=n_1+n_2+...+n_L$) is the concatenation of all $x_l$'s and, therefore, the entire search space of the problem. The objectives and constraints for each player are defined as $f^l: \real^{n}\mapsto\real$ and $g^l:\real^{n}\mapsto\real^{m_l}$. It is often assumed that no two players share any degrees of freedom (or decision variables) but we make no such assumptions in this work. To be precise, any player at level $l$ is a Stackelberg follower of all preceding players at level $1, 2, ...~ l-1$ and is simultaneously a Stackelberg leader for all players at level $l+1~...~L$. Like before, each player cares for their own objective and has their own constraints. To define the solution of the problem, we start with the concept of a rational reaction set for the final player L, $\phi^L(x_1, ... x_{L-1})$ defined as:

\begin{align*}
\phi^L(x_1, ... x_{L-1}) &:= \arg\min_{x_L} f^L(X) \\
&~~~~s.t.~~g^L(X) \ge 0
\end{align*}

Then, the rational reaction set for any player $l$, i.e. $\phi^l(x_1...x_{l-1})$ can be recursively defined as:

\begin{align*}
\phi^l(x_1...x_{l-1}) &:= \arg\min_{x_l, x_{l+1},~...~x_L} f^l(X) \\
&~~~~~~~~~s.t.~~g^l(X) \ge 0 \\
&~~~~~~~~~~(x_{l+1}...x_L) \in \phi^{l+1}(x_1...x_l)
\end{align*}

The solution to the entire problem is then:
\begin{align*}
\phi^1 &:= \arg\min_{x_1,~...~x_L} f^1(X) \\
&~~~~s.t.~~g^1(X) \ge 0 \\
&~~~~~~~~~~(x_{2}...x_L) \in \phi^2(x_1)
\end{align*}

It must be noted that in general, $\phi^l$ may not be a singleton, and therefore, there can be multiple local solutions for the problem.

These problems are not new and have been researched over the years in domains such as Economics, Optimal Control, Operations Research, and Decision Programming etc., for example, to model multi-stakeholder fund allocation, supply chain networks, inventory management, and power system security [Han et al., 2015, Yao et al., 2007, Fard and Hajiaghaei-Keshteli, 2018, Cassidy et al., 1971]. There are further generalizations of multilevel problems that include multiple players at each level (sometimes called a Multilevel Decentralized problem) who have equal deciding power amongst themselves but are, as before, followers for players above them and leaders for players below them. In this work, we restrict ourselves to multilevel optimization problems with only a single decision maker at each level and introduce a monte-carlo sampling based method to find solutions of such multilevel optimization problems. We also model a robust trajectory optimization problem and a generalized version of the toll-setting problem as multilevel problems and use our algorithm to find solutions for them. In summary, the main contributions of this paper are:

* A simple yet powerful monte-carlo method to solve multilevel problems.

* Modeling adversarial initial condition determination and nested toll-setting problem as multilevel optimization problems and obtaining their solutions via our proposed algorithm.

The remainder of this paper is structured as follows: In section <ref>, we explore some of the works related to such multilevel optimization problems, including some of the algorithms proposed to solve them. In section <ref>, we propose a stochastic algorithm to solve problems of this kind. Then, in section <ref>, we construct two such multilevel problems: a) a robust optimization problem of finding adversarial initial condition, and b) a nested toll-setting problem, and discuss the nature of their solutions. Then, in section <ref>, we apply this algorithm to solve a few problems from existing literature in addition to the constructed problems from section <ref> and compare the obtained solutions. In section <ref>, we perform empirical comparisons to study the convergence speed and computation time of the proposed algorithm. Finally, in section <ref> we pave the way for further research by outlining some of the possible improvements we envision in this domain and proceed to conclude the work with a brief recap.

§ LITERATURE REVIEW

Stackelberg games and Bilevel Optimizations are well-researched problems, and we refer readers to Dempe, 2020 in lieu of attempting a survey ourselves. Henceforth, we limit ourselves to works related to trilevel or general multilevel problems.

§.§ Linear Multilevel Problems
[Cassidy et al., 1971] first modeled the flow of resources between the federal, state, and municipal levels as a trilevel problem and provided a recursive dynamic algorithm for solving such problems. [Bard, 1984] later established stationarity conditions for trilevel linear optimization problems, generalized it to p-level stationarity problems, and devised a Cutting plane algorithm to solve them. [Ue-Pyng and Bialas, 1986] devised a hybrid method based on the K-th best algorithm and Parametric Complementary Pivot algorithm to solve trilevel linear problems. [Anandalingam, 1988] devised another method for solving trilevel linear problems by first obtaining and embedding the first-order necessary conditions (FONCs) of the third-level problem into the second-level problem, then obtained FONCs of thusly obtained problem and embedded it into the first-level problem. [Benson, 1989] investigated a specific case where linear multilevel problems are unconstrained and performed rigorous geometric analysis. Their major result was to show that the feasible solution set of such problems is a union of connected polyhedral regions. White, 1997 modified Bard, 1984 's method by changing the first step in their algorithm and claimed a qualitative improvement on the overall results.

§.§ Fuzzy Set / Goal Programming-Based Approaches

[Lai, 1996] considered a fuzzy set based algorithm to model and solve linear bilevel and multilevel problems. [Shih et al., 1996] later improved it to model problems that are not just hierarchical but also decentralized, or both, in nature. [Pramanik and Roy, 2007] modeled the multilevel problem as a fuzzy goal programming problem to solve it. [Zhang et al., 2010] presented a kth-best algorithm to solve linear trilevel programming problems and solved a constructed problem of annual budget allocation in a company with CEO, branch heads, and group supervisors.

§.§ Meta-heuristics based approaches

[Woldemariam and Kassa, 2015] developed a genetic algorithm based method to solve arbitrarily deep multilevel problems for bounded decision variables. [Han et al., 2016] devised a particle swarm optimization based method to solve bilevel problems and used it to solve a trilevel problem as well by embedding the stationarity conditions of the last level problem into the second level problem and converting the entire structure into a bilevel programming problem. At this point, we must also mention [Lu et al., 2016]'s survey of multilevel decision-making problems, which, although a bit dated, is an excellent resource for multilevel problems, algorithms, and applications developed until 2016.

§.§ Applications

[Han et al., 2017] used Vertex Enumeration method to solve a decentralized supply chain network involving manufacturers, logistic companies, and consumers modeled as a trilevel decentralized programming problem. [Fard and Hajiaghaei-Keshteli, 2018] modeled a multi-stakeholder supply chain problem as a trilevel problem and used five different meta-heuristic algorithms to solve them by solving each level in a turn-based fashion. They also later modeled a tire closed-loop supply chain network as a trilevel problem and solved it using a similar approach [Fard et al., 2018]. [Tilahun et al., 2012] developed a turn-based optimization strategy similar to [Fard and Hajiaghaei-Keshteli, 2018] to solve general Multilevel problems and later generalized it to solve fuzzy Multilevel, multi-objective problems with collaboration. [Tilahun, 2019]. [Tian et al., 2019] formulated a coordinated cyber-attack scenario as a trilevel problem and used the column and constraint generation method to obtain a solution. [Luo et al., 2020] modeled an Energy scheduling problem as a trilevel optimization problem and exploited its structure to obtain a closed analytical expression. [Laine et al., 2023] later developed a general algorithm to find solutions to Generalized Feedback Nash Equilibrium problems, which can be modeled as a Multilevel Stackelberg problem.

From the literature review, it is clear that multiple methods exist to solve trilevel problems, but only a few of these can be generalized to solve an arbitrarily deep multilevel problem. Even then, we find that each method has its own limitations. For instance, fuzzy set based methods ([Lai, 1996, Shih et al., 1996, Pramanik and Roy, 2007]) implicitly assume some degree of cooperation from lower levels, which is not an assumption that holds for every problem. Similarly, turn-based methods of [Tilahun et al., 2012, Tilahun, 2019, Fard and Hajiaghaei-Keshteli, 2018, Fard et al., 2018] are iterative best response algorithms that are more suited to find solutions to Nash equilibrium problems, and since they do not take into account the rational reactions of lower-level players, they do not converge towards the Stackelberg equilibrium. [Woldemariam and Kassa, 2015]'s genetic algorithm is quite promising, but it only works for bounded variables, which makes it inapplicable for a wide class of problems. Similarly, [Laine et al., 2023]'s algorithm is applicable only under assumptions of strong complementarity.

In light of these facts, we propose an algorithm in section <ref> that solves all of the outlined concerns above. Furthermore, we demonstrate in section <ref> that even though it's simple and intuitive, it outperforms the existing methods of similar nature. Compared to other algorithms, our proposed algorithm has the advantage that:

* It can handle problems with unbounded decision variables and, thus, is applicable to a wider class of problems.

* It can handle problems with non-differentiable objectives, as long as the final objective is differentiable.

* It can handle equality constraints present at the final level, unlike other Meta-heuristic algorithms, which fail to handle any equality constraints at all without any reformulations.

* It does not require any reformulations of the objective functions and, thus, can solve problems that can't be approached via KKT or Value function based reformulations.

* It's an anytime algorithm and can be tuned to obtain arbitrary accuracy at the expense of computation.

§ MONTE CARLO MULTILEVEL OPTIMIZATION (MCMO)

Some of the notations used in the algorithm are as follows:

\begin{align*}
L \in \mathbb{N} &: \text{Number of players}.\\
n_l \in \mathbb{N} &: \text{Number of variables for player }l\\
x_l \in \real^{n_l} &: \text{Variables that $l$-th player controls}\\
C^l &: \text{Feasible region for player }l\\
X \in \real^{n} &: \text{Concatenation of all $x_l$'s}\\
\end{align*}

\begin{align*}
C = \bigcap_{l=1}^{L} C^l &: \text{Feasible region for the problem }\\
x_s &: \text{Initially feasible point s.t. } x_s \in C \\
f^l &: \text{Objective function of player }l ~(\mathbb{R}^{n}\mapsto \mathbb{R}) \\
D^l := \real^{n_l} &: \text{Subspace spanned by}
\alpha^l \in \real^+ &: \text{Step size for player } l\\
N^l \in \mathbb{N} &: \text{Number of samples generated for player }l\\
M^l \in \mathbb{N} &: \text{Number of sampling iterations for } l
% \eta > 1 &: \text{Cooldown parameter}
\end{align*}

Apart from the notations above, we use some colloquial array notations as follows:

\begin{align*}
[~] &: \text{Empty array}\\
X[a:b] &: \text{Slice of X from index a to b inclusive}\\
X[a:end] &: \text{Slice of X from index a to the length}\\
&~~~ \text{of X inclusive} \\
X~.+b &: \text{A broadcasting summation operator.}\\
\end{align*}

MCMO is a sampling based algorithm. It iteratively refines any given approximate solution by generating samples in its neighborhood. These samples are successively passed down to each lower-level players, who generate samples of their own and pass them down to their lower-level players. This continues until the very last level, where a solver is used to obtain the solution for $x_L$ given the variables $x_1, ... x_{L-1}$ for the corresponding objective and constraints. Once these solutions are obtained, they are returned to upper-level players who evaluate them, select the best among them for their own objectives, and subsequently return them to their upper levels. At level 1, all returned solutions are evaluated, and the best among them is kept as the current estimate of the solution. In this way, MCMO acts as a gradient-free solver and does not require gradient information for any objective or constraint function except the last one. Similarly, since the last level is always solved by using a solver, MCMO can accommodate both equality / inequality constraints for that level so long as it's supported by the solver. MCMO is described in Algorithm <ref>. In essence, it takes an initially feasible point and continuously searches in its neighborhood for a better feasible and optimal point for a specified number of iterations. When the desired number of iterations is reached, MCMO returns a smoothed result from the last $k$ obtained iterates, as outlined in subsection <ref>.

$MCMO~(x_s, k)$
[1]
$X \gets x_s$
$P \gets [X]$
$i \in [1 ... maxiter]$
$X \gets OPTIMIZE(X, 1) ~or~ X$ # stick with same point if no better point found.
$P \gets P \cup X$
$\textbf{return} SMOOTHEN(P, k)$

The Optimize function defined in Algorithm <ref> takes as input an initially feasible point $x_s$ and a level $l$ ($=1$ for initial call). For the final player ($l=L$), this function uses a solver, IPopt [Wächter and Biegler, 2006] in this case, to optimize for the final objective $f^L$ subject to the constraints $C^L$. In all other cases, it generates $N^l+1$ random directions (including the zero direction) in the subspace $D^l$ to obtain new candidate points, which are then recursively passed to the optimizers of the lower-level player, i.e., $l+1$. These passed candidate points are then recursively perturbed by the lower-level players and returned. Out of all the returned values, player $l$ keeps the perturbed candidates that's best for its objective and satisfies its feasibility constraints. This process is repeated by the player $l$ for $M^l$ number of times, where, at the end of each such sampling iteration, it chooses the point that is the best among all obtained candidates in that iteration and uses it for the next iteration. At the end, it returns the final obtained best candidate to the upper player $l-1$. In the event that no feasible point can be found at any iteration, the last known best candidate point is retained and used for the next iteration. If no feasible point can be found even after $M^l$ iterations, the function returns $null$. In this way, this function can obtain solution for multilevel problems with arbitrary levels.

The algorithm uses three sub-procedures SOLVE_FULL, ARGMIN, and RAND_DIRECTIONS. These sub-procedures are intuitive, and thus, we only explain but do not explicitly outline them here. SOLVE_FULL takes, in order, an initial point, an objective, a set of constraints, and player level (to determine degrees of freedom to optimize on) and uses a solver to fully solve it to completion. Similarly, ARGMIN takes, in order, a list of candidate points, the player level $l$, and determines the best point according to the objective function $f^l$ ignoring any null points in the given list. Finally, RAND_DIRECTIONS generates $N^l$ random directions from a uniform hypercube (of length 1) centered at the origin in the subspace $D^l$.

$OPTIMIZE~(X, l)$
[1]
Require: $f^l, C^l, M^l, \alpha^l, N^l, D^l$
$X_R \gets null$
$l = L$
$X_R \gets SOLVE\_FULL(X, f^L, C^L, L)$
$X_R \not\in C^L$
return null
return $X_R$
$k \in \{1, ..., M^l\}$
# generate candidate points
$X_C \gets X ~{.+}~ \{ \alpha^l \cdot $
$~~~~~~~~~~~~~~~~~~~~~~~~~~~~~RAND\_DIRECTIONS(N^l, D^l) \cup \textbf{0} \}$
$x \in X_C$
$x \gets OPTIMIZE(x, l+1)$
$x \not\in C^l$
$X_R \gets ARGMIN(X_R, x, l)$
$X \gets X_R ~or~ X$
return $X_R$

Ideally, MCMO should be used with a high number of samples and sampling iterations, i.e. $N^l, M^l$ to obtain accurate results, as only by doing so can we solve all lower levels to completion before optimizing any upper-level problem. But this can result in a lot of computational overhead, as outlined in section <ref>. So, in practice, we select a reasonable $N^l, M^l$ for each level while still keeping the problem computationally tractable. But this will result in stochastic estimates (as opposed to true solutions), which is precisely what limits MCMO to an approximate algorithm.

§.§ Initialization

MCMO requires that an initially feasible (not necessarily optimal) point $x_s \in C$ be provided. A viable option to achieve such an initially feasible point is to solve the following problem:

\begin{align*}
x_s = &\arg\min_{X}~~ 0\\
&~~~~~~~~ s.t. ~X \in C
\end{align*}

However, a heuristic to achieve a reasonably optimal starting point for a non-trivial problem is to take the weighted sum of the objectives. Which is to say, we solve the following optimization problem to obtain such an initially feasible point:

\begin{align*}
x_s = &\arg\min_{X}~~ \sum_{l=1}^{L} w_l f^l(X)\\
&~~~~~~~~~~~ s.t. ~X \in C
\end{align*}

Where, $w_l$'s are chosen as required. This heuristic yields better starting points in cases where the true solution lies closer to the pareto front of the involved objective functions.

§.§ Smoothing

Since MCMO is a stochastic algorithm, it can only provide approximate solutions as all the lower-level problems are not completely solved. This is especially true when the number of samples generated ($N^l$) or the number of sampling iterations ($M^l$) are too low and $\alpha^l$ is high. Therefore, at the end of the algorithm, the last $k<maxiter$ points are used to obtain a more stable approximation of the equilibrium point by using a smoothing scheme. The choice of smoothing scheme may depend upon the problem, but in this work, we use the following scheme:

Best Objective Smoothing Scheme: $X^*$ is approximated as $X^*=\arg\min_{X \in \{X_1, X_2, ... X_k\}} f^1 (X)$. Where $f^1$ is the objective function of the first player. This scheme is guaranteed to produce a feasible point (since all $X_1, X_2, ... X_k$ are feasible, as shown in subsection <ref>).

§.§ Practical Considerations

The performance of the algorithm is reliant on the number of samples $N^l$ generated per level, number of sampling iterations $M^l$, choice of $\alpha^l$, and the number of iterations $maxiter$. In general, more samples and sampling iterations would improve the accuracy of the solution, but at the expense of computation costs. Similarly, a large $\alpha^l$ may prevent convergence, whereas a low $\alpha^l$ would delay it. An appropriate way of running MCMO is thus to start off with low $N^l, M^l$ and high $\alpha^l$ and then fine-tune the result with a lower value of $\alpha^l$ and a higher sample size $N^l$ and iteration $M^l$ to the desired accuracy. Due to the nature of multilevel optimization problems, lower-level players must be provided with greater deciding powers than any upper-level players. This is especially true when the degrees of freedom are shared between the upper and lower level players. Consider for example, the problem

max_x  x

s.t.  x ∈      min_x x

s.t. x ∈[l, u]

In this case, the solution for this problem for $x\in[l, u]$ is $l$. In terms of games, when the first-level player chooses any $x=x'$, the second-level player will choose $x=l$, overriding any choice of the variable $x$ made by the upper-level player.
Therefore, to achieve true solutions, $\alpha^l,$ $N^l, M^l$ for each subsequent level should be increased. Additionally, the choice of $N^l$ should also consider the degrees of freedom. If a player has control of two variables, they must be allowed to sample more directions than if they only had one decision variable. This ensures that the sampling is fair for all levels.

However, for simple problems, it may also be desirable to use the same $\alpha$ per player for maintaining a lean parameter space. And if bounds on the player's variables are known, it can guide the choice of $\alpha$.

§.§ Computation Time

MCMO is a recursive sampling based algorithm and thus, its computation time increases exponentially with each additional level. Furthermore, the computation time will also depend upon the parameters $N^l, M^l, maxiter$ and the nature of the problem itself. While parallelizing the implementation may provide speedups, for this work, we do not attempt such efforts and have left it for future improvements. A detailed empirical analysis of computation time can be found in section <ref>.

An implementation of the algorithm can be found on GitHub (https://github.com/VAMPIR-Lab/MCMO).

§.§ Proofs

This section presents proofs for the feasibility and convergence of the MCMO algorithm.

§.§.§ Proof of Feasibility

Any non-null point X returned from a function call of the form $X=~$OPTIMIZE($\cdot$ , l) is feasible for level l, i.e., $X \in C^l$.

For the final level $l=L$, this is easy to see from Algorithm <ref> lines 4–8. If a point is infeasible for that level, the if condition on line 4 causes a null return. Otherwise, a feasible point is returned in line 7. For levels $l\neq L$, a non-null result can only be returned if $X_R$, which is initially null, is set with a non-null value $x$ in line 18. But line 18 can only execute if the feasibility condition of line 15 was satisfied, which means that the returned non-null value $X\in C^l$.

Any non-null point X returned from a function call of the form $X=$ $OPTIMIZE(\cdot, l)$ is always obtained from a lower-level function call of the form $OPTIMIZE(\cdot, l+1)$ for $l < L$.

Since $l \neq L$, following arguments similar to lemma <ref>, it must have been set by line 18. But any such point is clearly obtained in line 14 by function call of the form $OPTIMIZE(\cdot, l+1)$. Hence, this is true.

We can now prove the following claim:

Each iteration in MCMO function obtains a feasible point.

From lemma <ref>, we know that any non-null point obtained from function of the form $OPTIMIZE(\cdot, 1)$ is obtained from $OPTIMIZE(\cdot, 2)$, $OPTIMIZE(\cdot, 3)$, and so on until $OPTIMIZE(\cdot, L)$. Similarly, we also know from lemma <ref> that any non-null point thus obtained must be feasible for levels 1, 2, ... $L$-1, and L. Therefore, any non-null point obtained from an iteration of the MCMO algorithm is feasible for all levels. Furthermore, if a null point is obtained at any point, MCMO retains the last non-null point, which is either $x_s$, an initially feasible point, or another non-null point previously obtained in iteration that has already been shown to be feasible.

§.§.§ Proof of Convergence

Any analytical reasoning for general multilevel problem is decidedly hard, and for stochastic or meta-heuristic algorithms, the difficulty only increases. Thus, we only present an asymptotic proof of convergence for a narrow class of problems that satisfy the following simplifying assumptions:

* The rational reaction set $\phi^l(x_1, ..., x_{l-1})$ (as defined in section <ref>) for player $l$ is a point-to-point map, i.e., all rational reactions are unique for given upper-level decisions.

* A solution exists for the given problem, and the solver used for the final level can always find solutions when they exist.

In general, assumption 1 may not be valid but may hold if the upper-level constraints are restrictive enough or if the topmost objective is strongly convex and we want to solve an optimistic multilevel optimization problem, i.e. lower levels cooperate with the topmost player for ambiguous rational reactions. Furthermore, this is a simplification that multiple analytic treatments of this problem [Liu, 1998, Woldemariam and Kassa, 2015] have made as arguing about the problem in general is intractable.

Under our assumption, for any multilevel Stackeblerg problem, the optimization that player $l$ solves, say $P^l(x_1, ... x_{l-1})$, condenses to:

\begin{align*}
P^l(x_1, ... x_{l-1}) := \min_{x_l} &f^l(x_1, ..., x_l, \phi^{l+1}(x_1, ..., x_l)) \\
&~~~~s.t.~~g^l(X) \ge 0 \\
\end{align*}

$OPTIMIZE(\cdot~, L)$ solves $P^L(x_1, ... x_{L-1})$.

For the last level, i.e., $l=L$, this function uses a solver to obtain the solution. Since it's assumed that a solution exists and that the solver can find it, this is trivially true.

If $OPTIMIZE(\cdot~, l+1)$ solves $P^{l+1}(x_1, ... x_l)$, then $OPTIMIZE(\cdot~, l)$ solves $P^l(x_1, ... x_{l-1})$ given $N^l, M^l \to \infty$

Since infinite samples are assumed with infinite sampling iterations and it's also assumed that a unique reaction (say $x_l^*$) exists for $P^l$, we claim that the sampling process would eventually converge towards $x_l^*$. To show that this is indeed true, we first assume that the sampling does not converge towards the optimum $x_l^*$. This can only mean one of the following:

* The algorithm cycles between points $x_l^1, x_l^2, ... x_l^i$. But this must mean that $v(x_l^1) > v(x_l^2) > ... > v(x_l^i) > v(x_l^1)$, which is a contradiction. Here, we define $v(x_l):= f^l(x_1, ..., x_l, \phi^{l+1}(x_1, ... x_l))$.

* The algorithm gets stuck on some $x_l$ and no $x_l'$ exists in its neighborhood such that $v(x_l')<v(x_l)$ and $x_l'$ satisfies appropriate constraints. However, this, by definition, is a local optimum for the player $l$ and thus, by our assumption, is the same as $x_l^*$ and results in contradiction.

MCMO eventually converges upon the unique solution.

Under our framework, the overall problem reduces to $P^1$. From lemma <ref> and <ref>, we have a proof by induction that MCMO solves $P^1$ when $\forall l, N^l, M^l \to \infty$ by calling $OPTIMIZE(\cdot~, 1)$

§ SOME MULTILEVEL PROBLEMS

§.§ Adversarial Initial Condition (AIC) determination problem

We can loosely define a Trajectory as a continuous path in some space. In robotics and control, such paths are generally produced from some initial conditions (start point, environment, etc.) by a set of rules or functions, usually called a policy.
This problem is related to finding a worst-case initial condition for any given policy. The worst-case being an initial point from where, if a trajectory is generated according to such a policy, it ends up a) bringing the trajectory as close to touching the obstacle as possible, and b) increasing the length cost of the trajectory. Figure <ref> depicts the problem we construct here.

We consider a 2D plane to be our environment. The blue circular region is the feasible region $\chi \subset \real^2$ where any start point $x \in \real^2$ is allowed to reside. A fixed and known policy $\Pi$ then generates a trajectory $\Tau = \Pi(x) = \tau^0, \tau^1, ...~, \tau^i \in \real^2, ~ \tau^0 = x$ up to the finishing line $D \in \real$ using the start point such that some cost $f(\Tau) \in \real$ (modeled here as the horizontal length of the trajectory i.e. $f(\Tau) = D-\tau^0_1$) is minimized and certain feasibility conditions for each trajectory points are satisfied i.e. $g(\tau^i) \ge 0~\forall \tau^i \in \Tau$.

In this example, the condition of feasibility for a trajectory $\Tau$ is that all trajectory points be outside the obstacle region $\mathcal{O} \subset \real^2$. Modeling $\mathcal{O}$ as a circle centered at $o$ with radius $r$, our feasibility condition for each trajectory point becomes: $g(\tau^i) = ||o-\tau^i||^2 - r^2 \ge 0 ~~ \forall \tau^i \in \Tau$. The problem that we consider in this work is to find an adversarial initial point $x^a$ such that for any given policy $\Pi$, the generated trajectory $\Tau^a = \Pi(x^a)$ is as close to infeasibility and sub-optimality as possible. The rationale being that, with such obtained point, we could iterate our policy to improve it under even the most adverse initial conditions. We do not attempt policy training in this text and have left it for future work.

The environment for the AIC problem. $\Tau$ is a sampled sinusoidal trajectory.

We can model this problem as a trilevel game as follows:

\begin{align}
\label{eq:trilevel}
&\max_{x,T\in\real}~~f(\Tau) \\
&~~~~ s.t~~x,T\in\arg\min_{x, T}~ T\\
&~~~~~~~~~~~~~~~~~~~ s.t.~~ x \in \chi \\
&~~~~~~~~~~~~~~~~~~T \in \arg\max_{T} ~ T\\
&~~~~~~~~~~~~~~~~~~~~~~~~~~~~~s.t~T \ge 0\\
&~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ \Tau = \Pi(x)\\
&~~~~~~~~~~~~~~~~~~~~~~~~~~g(\tau^i) \ge T;~~ \forall \tau^i \in \Tau
\end{align}

An interpretation of the trilevel problem is as follows: The first player wants the initial point $x\in \chi$ for the trajectory to maximize the cost of the trajectory $f(\Tau)$. The second player, whereas, wants to bring $T$ close to 0 by manipulating $x$. But, $T$ is the minimum of all feasibility scores $g(\tau^i)$ i.e., the point closest to violation. So when $T\to 0$, the closest trajectory point to the obstacle becomes even closer, and as a result, the trajectory $\Tau$ touches the obstacle $\mathcal{O}$. Here, the first and second players share the same degree of freedom, i.e., the variable $x$. Generally, in multilevel games such as these, non-overlapping d.o.f's are considered. But, as we mentioned previously, we make no such assumptions and design a general algorithm that can handle all such scenarios.

§.§ Nested Toll-Setting problem

A toll-setting problem is a well-known bilevel optimization problem where a toll-setter decides a toll amount for a road segment. Since they want to maximize the total income, they can neither set the toll too high, or drivers will avoid the road segment due to exorbitant fees, nor set the toll too low, or their total income will decrease. We refer readers to [Labbé et al., 1998] for a more detailed treatment of this problem. Instead, we focus on a generalization of this problem, i.e., the Nested Toll-Setting problem, as shown in figure <ref>.

The nested toll-setting problem. The numbers in black above a road segment represent the percentage of traffic on that segment. The numbers in red below a road segment represent the cost of taking that segment. Costs are sum of toll price, congestion, and other extra factors.

We consider two toll stations $T_1$ and $T_2$ established to oversee their respective tolled segments (red for $T_1$, and yellow for $T_2$). Any vehicle arriving at $T_1$ has the option to either take the tolled segment (red) by paying $t_1$ cost per unit traffic or take the non-tolled segment (black) for free. We assume that $p_1$ percentage of the original fleet takes the red tolled segment. Similarly, any vehicle arriving at $T_2$ has the option to either take the tolled segment (yellow) by paying $t_2$ per unit traffic or take the non-tolled segment. We assume that $p_2$ percentage of the original fleet takes the yellow tolled segment and $p_3$ percentage of original fleet takes the final free segment. It must be clarified that $p_1, p_2,$ and $p_3$ represent the percentage of fleet that first arrives at $T_1$ establishing $p_1 + p_2 + p_3 = 1$.

From the perspective of the fleet, the cost of travelling through any segment is the sum of a) the toll on the segment, b) congestion on the segment, and c) additional costs associated with the segment. For the purpose of this problem, we establish the congestion cost for any segment to be $\sigma \cdot p$, where $p$ is the percentage traffic on the segment and $\sigma$ is a constant. This is to say that the congestion cost increases linearly with the traffic on the segment. We further simplify this problem by setting $\sigma=1$, thereby setting the congestion cost at each segment equal to the traffic percentage at that segment. Finally, we assume that none of the road segments have any additional costs except for the final free segment, which has an extra cost $D$. This extra cost could theoretically model road length, road conditions, traffic lights, or a myriad of other factors. For this problem, we allow $D$ to be less than 0, allowing it to model a reward or a subsidy as well. Figure <ref> shows the costs associated with each road segment below the segments in red.

Once the toll has been set for any segment, the fleet decides to divide a certain percentage of its traffic to the tolled road or the free road to minimize its total cost across each road segment. Now if, at each toll station, the fleet makes a greedy decision, i.e., deciding whether or not to take the tolled road without considering any future toll stations on the way, then this problem can be written as the following trilevel problem:

max_t_1, t_2,p_1,p_2,p_3     p_1 ·t_1 + p_2 ·t_2

s.t. t_1, t_2    ≥0

p_1,p_2,p_3   ∈min_p_1,p_2,p_3 p_1 ·(p_1 + t_1)

+ (p_2 + p_3)^2

s.t. p_1 ∈[0, 1]

p_2,p_3∈min_p_2, p_3   p_2 ·(p_2+t_2)

+ p_3 ·(p_3 + D)

s.t. p_2 ∈[0,1]

p_3 ∈[0,1]

p_1+p_2+p_3 = 1

Here, the first level corresponds to the toll-setter who decides on $t_1, t_2$ to maximize their total income. The second level is the fleet's decision at station $T_1$. It chooses the percentage of traffic to balance total congestion costs and toll costs. The third level is the remaining fleet's decision at station $T_2$ for the same.

§.§.§ Nature of the solution

It can be argued with relative ease that for a very high $D$, it's beneficial for the toll-setter to redirect all traffic to $T_2$ whereas for a very low $D$, the toll-setter is better off exacting all tolls from $T_1$ instead. In fact, there are two known equilibrium points $X = (t_1, t_2, p_1, p_2, p_3)$ for this problem for different values of $D$ (see Appendix <ref>):

* At $D = \olsi D = 6, \olsi X = (\ge 2, 4, 0, 1, 0)$

* At $D = \underline D = -1.5, \underline X = (1, \ge 0, 0.25, 0, 0.75)$

§ EXPERIMENTS AND RESULTS

In this section we solve some existing multilevel problems from the literature in addition to our constructed problems, i.e., the adversarial initial condition (AIC) problem and the nested toll-setting problem outlined in section <ref> using MCMO algorithm of section <ref>. To keep our parameter space restricted and the experiments simple, we run each of these problems with the same value of $\alpha$ for different levels. Furthermore, we also set all $M^l=1$ and instead setup our algorithm based solely on $\alpha^l$ and the number of samples $N^l$. All examples have been run on a personal computer with an Intel Core i5 8400 processor with a 2.8GHz frequency and 32 GBs of DDR4 RAM.

§.§ Solving AIC using MCMO

We now solve the AIC problem for two policies, as described below.

§.§.§ Linear Policy: $\Pi^l$

We define the linear policy $\Pi^l: \real^2 \mapsto \real^2$ as follows:

$$\Pi^l([x_1, x_2]^T) = [x_1+\delta, x_2]^T$$

Where $\delta \in \real$ is a step-size. Intuitively, this policy takes a point $x^i$ and generates a point $x^{i+1}$ by stepping $\delta$ distance in the $x_1$ axis while leaving $x_2$ unchanged, i.e., a horizontal line parallel to the $x_1$ axis.

§.§.§ Non-Linear Policy: $\Pi^n$

We define $\Pi^n : \real^2 \mapsto \real^2$ as follows:

\begin{align*}
\Pi^n([x_1, x_2]^T) = [ &x_1+\delta, \\&x_2 + A \left( \sin( B (x_1 + \delta)) - \sin( B (x_1)) \right)]^T
\end{align*}

Where $\delta \in \real$ is a step-size, $A \in \real$ is an amplitude parameter, and $B \in \real$ is frequency parameter. This generates a sinusoidal trajectory parallel to the $x_1$ axis.

§.§.§ Setup for AIC

For both of the trajectories, we apply MCMO to obtain adversarial points for different placements of the obstacle circles of radius r = 2. For all experiments, our feasible region is a circle centered at $[5, 5]^T$ with a radius of 5 r. The number of trajectory points is fixed at $N^\tau = 20$, and the destination plane is set to $D = 20$. For both policies, step-size $\delta$ is set to 1, and for non-linear trajectory $\Pi^n$, $A, B$ are set to $0.5, 3$, respectively. We apply MCMO for a maximum of 150 and use the best objective smoothing scheme with 10 final samples. Similarly, the step parameter "alpha is set to 3. For both policies, $N^1$ was chosen to be $2$ and $N^2$ was chosen to be $10$. In general, $N^2>N^1$ is in accordance with subsection (<ref>), but in addition, player 2 has more degrees of freedom as compared to player 1, and furthermore, both player 1 and player 2 share two degrees of freedom $(x_1, x_2)$, so no matter what player 1 chooses, it is modified by player 2, so player 1 has very little influence to begin with. As discussed previously, we set $M^1=M^2=1$. While initializing, we found that the weights $w_1 = 10^5, w_2 = 10^{-5}, w_3 = 1$ (<ref>) gave us feasible starts. In general, this will always depend on the problem being solved. The initial points produced for linear policy are shown in figure <ref>, and those for nonlinear policy are shown in figure <ref>.

Adversarial Initial points for linear trajectories for obstacle center $o=(15, 13)$ (top-left), $o=(12, 9)$ (top-right), $o=(15, 5)$ (bottom-left), and $o=(12, -3)$ (bottom-right). Time taken for the solutions is, in order, 220.8 seconds, 269.47 seconds, 262.9 seconds, and 267.39 seconds. Differences in timing indicate that the final level solver converged quickly for some instances of the problem. Ideal points are as left as possible and are either touching the obstacle or come as close to touching it as possible. Red lines indicate the path of solution in the feasible region.

Adversarial Initial points for nonlinear (sinusoidal) trajectories for the obstacle center are $o=(15, 13)$ (top-left), $o=(12, 9)$ (top-right), $o=(15, 5)$ (bottom-left), and $o=(12, -3)$ (bottom-right). The time taken for the solutions, in order, is 2272.13 seconds, 2624.19 seconds, 2408.66 seconds, and 2356.23 seconds. Ideal points are as left as possible and are either touching the obstacle or come as close to touching it as possible. Red lines indicate the path of the solution in the feasible region. All the solutions obtained are quite close to optimality. The top-left instance may not look optimal but the phase of the sinusoid and our sampling strategy may go counter to our intuition.

§.§ Discussion on Results

As can be seen, in all cases, MCMO generates proper adversarial initial conditions for this problem. For linear policy (figure <ref>), except the bottom-right setup, all other instances of the problem achieve optimal results. For the bottom-right instance, although the obtained point is not optimal, the error is $\approx 7.5\%$ which is not at all unreasonable considering the stochatic nature of the algorithm. However, accuracy can be further increased to desired bounds by running MCMO with a higher number of samples and sampling iterations and lower values for $\alpha$ for further iterations. The path taken by the solution at each iteration is traced by the red line. Unsurprisingly, for problems where the initial feasible solutions were closer to optimality, the algorithm converged to the answer in very few iterations. For problems where the initially feasible solution was not close to optimality, the path appears repetitive and chaotic, eventually converging to the answer, but it must be taken into account that the plotted path is a projection $(x_1, x_2)$ of the true decision space $(x_1, x_2, T)$.

For non-linear policy (figure <ref>), almost all of the instances converge to the optimum. While the true trajectory (represented by the sinusoid) does intersect the obstacle, it is to be expected because our problem formulation has been for the discrete samples to begin with, which incidentally behave as expected. Furthermore, it may also appear that the top-right instance of the problem is not optimal, as the generated trajectory is not as close to the obstacle as possible. But owing to the fact that the sinusoid's phase depends upon the initial point and taking into account our sampling strategy, moving the point to the top does not, in fact, bring the trajectory any closer to the obstacle.

§.§ Solving Nested Toll-Setting Problem using MCMO

A comparison of the known solutions for the edge cases with the solutions obtained by MCMO is tabulated in Table <ref>. The parameters used for both of the instances of the problem are $N^1=7, N^2=7, X_s=[0, 0, 1, 0, 0], \alpha=0.15, maxiter=100$. Smoothing scheme used is best objective smoothing with $k=10$. From the table, we can see that MCMO achieves results with error (w.r.t $f_1$) of $0.025\%~ \text{and} ~1.76\%$ respectively for parameters $D=6,-1.5$ taking, respectively, 119.11 and 129.73 seconds. The achieved results are quite satisfactory but can be made more accurate by decreasing the step sizes $\alpha$ and increasing the number of samples and sampling iterations as required.

$D$     $X^*$     $t_1$     $t_2$     $p_1$     $p_2$     $p_3$     $f_1$
$6$     $\olsi X$     $\ge$ 2     4     0     1     0     4

$X_{M}$     2.387     4.032     0.006    0.989     0.005     4.001
$-1.5$     $\underline X$     1     $\ge$0     0.25     0     0.75     0.25

$X_{M}$     1.281     0.15     0.199     0     0.801     0.254

A comparison of solution obtained via MCMO $X_M$ with known analytical solution $X^*$ for the Nested Toll-Setting Problem for different parameters $D=\olsi D = 6, \text{and} ~D = \underline D = -1.5$. The achieved errors are $0.025\%$ and 1.76% respectively.

§.§ Numerical Examples from the Literature

The following problem derived from [Sinha, 2003] is a trilevel linear problem defined as:

\begin{align*}
\label{eq:sinha}
&\max_{x_1, x_2,x_3,x_4}~~7x_1 + 3x_2 -4x_3+2x_4 \\
&~~~~ s.t.~x_3,x_4\in\arg\max_{x_3, x_4}~ x_2 + 3x_3 + 4x_4\\
&~~~~~~~~~~~s.t.~x_4\in\arg\max_{x_4} ~ 2x_1+x_2+x_3+x_4\\
&~~~~~~~~~~~~~s.t.~x_1+x_2+x_3+x_4\le 5 \\
&~~~~~~~~~~~~~~~~~~x_1+x_2-x_3-x_4\le2 \\
&~~~~~~~~~~~~~~~~~~x_1+x_2+x_3 \ge 1\\
&~~~~~~~~~~~~~~~~~~x_1-x_2+x_3+2x_4\le 4\\
&~~~~~~~~~~~~~~~~~~x_1,x_2,x_3,x_4 \ge 0
\end{align*}

The optimum $f_1^* = 16.25$ for this problem is reported at $(2.25, 0, 0, 0.25)$. MCMO obtains a result of $16.145$ at $(2.205, 0.06, 0, 0.265)$ when run with the following parameters: $N^1=6, N^2=3,M^1=M^2=1, X_s=[0.4, 0.4, 0.4, 0.4], \alpha=1, maxiter=100$. The sample size was chosen owing to the difference in the number of variables, while the feasible set was deduced by observation. The smoothing scheme used is the best objective smoothing with $k=10$. The relative error in objective values for this example is $< 1\%$. The time taken to obtain the solution is 59.26 seconds.

The second problem is taken from [Tilahun et al., 2012] and is defined as:

\begin{align*}
&\min_{x,y,z}~~-x + 4y \\
&~~~~ s.t. ~x+y \le 1\\
&~~~~ y,z\in\arg\min_{y,z}~ 2y+z\\
&~~~~~~~~~~s.t.~ -2x+y \le -z\\
&~~~~~~~~~~~z\in\arg\min_{z} ~ -z^2+y\\
&~~~~~~~~~~~~~~~~~~s.t.~~z\le x; x\in[0, 0.5]; y\in[0,1]; \\
\end{align*}

The reported optimum $f_1^*=-0.5$ is at $(0.5, 0, 0.0095)$ whereas MCMO obtains $f_1=-0.498$ at $(0.498, 0, 0.498)$ when run with the following parameters: $N^1=5, N^2=5, M^1=M^2=1, X_s=[0, 0, 0], \alpha=0.2, maxiter=100$. The choice of $\alpha$ was guided by the bounds on the decision variables, and the initial feasible solution was obtained via observation. The smoothing scheme used is the best objective smoothing with $k=10$.
The obtained minimizer disagrees with the reported minimizer, but it can be seen that the reported minimizer, i.e., $(0.5, 0, 0.0095)$ is incorrect as opposed to the actual minimizer, i.e., $(0.5, 0, 0.5)$ because once $x,y=0.5, 0$ are chosen, $z$ can be clearly increased (upto $x$) by the last player to achieve further minimality. Moreover, [Woldemariam and Kassa, 2015] agrees with our results on the same problem and reports $f_1=-0.4929$ at $(0.4994, 0.0016, 0.4988)$. For this problem, MCMO achieves a relative error of $< 1\%$ in 20.85 seconds.

§ COMPARISONS

§.§ Comparisons with Existing Works

We compare the results obtained in subsection <ref> for the Nested Toll-Setting problem with some of the existing methods from the literature. We chose [Tilahun et al., 2012] and [Woldemariam and Kassa, 2015] as baselines because these methods have been proposed for arbitrarily deep multilevel optimization problems as well. Since none of these methods are capable of solving problems with equality constraints, we have to reformulate the Nested Toll-Setting problem to remove the equality constraint as follows:

\begin{align*}
% \label{eq:newtoll}
&\max_{t_1, t_2,p_1,p_2}~~p_1 \cdot t_1 + p_2 \cdot t_2 \\
&~~~~ s.t.~ t_1, t_2 \in [0,10] \\
&~~~~ p_1,p_2 \in\arg\min_{p_1,p_2}~~p_1 \cdot (p_1 + t_1) + (1-p_1)^2\\
&~~~~~~~~~~~~~~~~~~ s.t. ~p_1 \in [0, 1] \\
&~~~~~~~~~~~~~~~~~~~~~~~~p_2\in\arg\min_{p_2} ~ p_2 \cdot (p_2+t_2) ~+ \\
&~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~(1-p_1-p_2) \cdot (1-p_1-p_2 + D)\\
&~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~s.t.~p_2 \in [0,1-p_1]\\
\end{align*}

We note that, in general, it may not always be possible to reformulate a problem to remove a particularly tricky equality constraint. For both methods, we generate an extremely high number of samples per iteration, i.e., $10^6$ for the entire decision space, and run both algorithms for 100 iterations. For $D=-1.5$ [Tilahun et al., 2012]'s method obtains $f_1=2.53\times 10^{-7}$ at $X=(27.65, 36.93, 1.21\times10^{-9}, 5.95\times10^{-9})$ and for $D=6$, it obtains $f^1=2.43\times10^{-6}$ while taking 875.2 seconds and 883.41 seconds respectively. These results have a very high error compared to the theoretical best, and the reason for this is that this algorithm does not truly solve a Multilevel Stackelberg problem at all. It's instead an Iterative Best Response type algorithm, which only works when finding the Nash Equilibrium of a problem.

[Woldemariam and Kassa, 2015]'s approach only works for bounded decision variables, so we add two additional constraints, i.e., $t_1\in[0,10], t_2\in[0,10]$. For $D=6$, it obtains $f_1=4.75$ at $X=(9.64, 5.13, 0.066, 0.8)$ in 500.76 seconds, which overestimates the theoretical maximum by the relative error of $18.75\%$ and for $D=-1.5$, it obtains $f_1=0.729$ at $X=(9.8, 6.72, 0.066, 0.0113)$ which has a relative error of $191.6\%$ in about 465.02 seconds. Even though this method works much better than the former, it still ends up overestimating the leader's objective most of the time. This is because of the update rule of this algorithm, which only ever changes the obtained solution if it's better than the previous one for just the leader. So if a solution with high complementary error but better leader objective is acquired, it's always kept regardless of whatever may be found in subsequent iterations.

§.§ Timing Comparison for arbitary levels

To compare the time required by MCMO to solve any given problem against its complexity, we introduce the following arbitrarily multilevel problem parameterized by $w\in(\real^+)^{n}$.

\begin{align*}
% \label{eq:nridgelin}
&\min_{x_1,...,x_n\in\real^n}~~ \norm{\column{x_1 \\ x_2 \\ \vdots \\x_n}-\column{w_1\\ w_2\\ \vdots\\ w_n}}\\
&~~~~~~~~~~~~ x_2, ...x_n \in \arg\min_{x_2,...,x_n\in\real^{n-1}}~~\norm {\column{x_2 \\ \vdots \\x_n} - \column{w_2\\ \vdots\\ w_n}}\\
&~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ s.t.~~x_2 \le x_1\\
&~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ \vdots\\
&~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~x_n \in\arg\min_{x_n\in\real} ~ (x_n-w_n)^2\\
&~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~s.t. ~~ x_n\le x_{n-1}
\end{align*}

While this problem generally degenerates to a single-level problem of the form $min_x \norm{x-w}; x_i\le x_{i-1} \forall i$, it's still an ideal problem to test our algorithm for the fact that the dimension of the decision variable $x$ increases linearly with the level. When solved with the parameters $M^i=1, \alpha_i=0.25 \forall i$, for 50 iterations each for different values of sample sizes per level $N^i$, the obtained timing results are tabulated in Table <ref> and the corresponding graph is shown in Figure 5. The execution time grows exponentially as the levels increase, which is as expected of a recursive algorithm.

4|c|Time (s)
Levels    N=3   N=4   N=5   N=6
2    1.33    1.68    1.71   2.33
3    3.84    5.52    7.71   11.86
4    10.03    21.57    39.41   77.194
5    29.84    101.81    433.91   1356.82
6   97.15    1230.33    7,463.52   10938.59

Time taken as Levels increase for different sample sizes $N$

title=Problem Complexity vs Time for different $N^i$,
xlabel=Number of Levels,
ylabel=Time for 50 Iterations [s],
xmin=0, xmax=7,
ymin=0, ymax=5000,
ytick=0,1000, 2000, 3000, 4000,5000,6000,7000,8000,9000,
legend pos=north west,
grid style=dashed,
(2, 1.33)(3, 3.84)(4, 10.03)(5, 29.84)(6, 97.15)
(2, 1.68)(3, 5.52)(4, 21.57)(5, 101.81)(6, 1230.33)
(2, 1.71)(3, 7.71)(4, 39.41)(5, 433.91)(6, 7463.52)
(2, 2.33)(3, 11.86)(4, 77.194)(5, 1356.82)(6, 10938.59)
$N^i=3$, $N^i=4$, $N^i=5$, $N^i=6$

Execution time as problem complexity increases for different sample sizes. The obtained graph is exponential as the problem increases linearly, which is as expected.

§.§ Accuracy Comparison

In general, we expect the accuracy of MCMO to increase as the number of samples per level $N^i$ increases. For this experiment, we use a five level version of the problem introduced in the previous subsection with a randomly generated $w = (3, 8, 7, 7, 3)$. We fix $M^i=1, \alpha^i=0.25$ and start with the initial guess $x_s=(0,0,0,0,0)$ and plot the convergence of the algorithm per iteration in Figure 6. As can be noticed, as the sample sizes increase, the convergence of the algorithm does increase, but it gets capped beyond a certain point because of the step size $\alpha^i$. The effect of increasing the step size for the same problem by keeping the sample sizes fixed at $N^i=6$ for different $\alpha^i=0.25, 0.5, 1$ is shown in Figure 7. It can be observed that increasing $\alpha^i$ for an appropriate number of samples increases the convergence speed to the optimum. Once the optimum is approached the convergence plateaus. This demonstrates that MCMO is stable, i.e., once it approaches the neighborhood of a stable solution, it remains there (provided enough samples are taken at each level).

title=Convergence for different $N^i$,
ylabel=Objective of first player $f^1$,
xmin=0, xmax=25,
ymin=70, ymax=200,
legend pos=south west,

$N^i=2$, $N^i=3$, $N^i=4$, $N^i=5$, $N^i=6$

Convergence for different sample sizes. After a certain threshold, step size $\alpha^i$ caps the performance.

title=Convergence for different $\alpha^i$,
ylabel=Objective of first player $f^1$,
xmin=0, xmax=25,
ymin=0, ymax=200,
legend pos=north east,

$\alpha^i=0.25$, $\alpha^i=0.5$, $\alpha^i=1$, True

Convergence for different step sizes for $w=(3,8,7,7,3)$. It can be shown that the minimum for this problem is at $X=(6.25, 6.25, 6.25, 6.25, 3)$ with leader's objective $f^1=14.75$ (shown in the plot by the cyan horizontal line)

§ CONCLUSION AND FUTURE WORK

Stackelberg games arise in many real-world scenarios, and conversely, many interesting economic, control, and other causal phenomena can be naturally modeled as Stackelberg games. Multilevel Stackelberg games provide a further generalization that expands the perimeter of interesting interactions that can be modeled by such rules. However, the difficulty involved in solving them is non-trivial and can present a major challenge to those who seek to model and solve these kinds of problems. In this paper, we introduced two such example problems that can effortlessly be modelled using multilevel formulation, i.e., a) the Adversarial Initial Condition determination problem, where we find a challenging initial condition for any provided policy, and b) the Nested Toll-Setting problem, which is a generalization of the famous Bilevel Toll-Setting problem. We then presented MCMO, a stochastic algorithm that can be used to solve problems involving an arbitrary number of leaders and followers (i.e., arbitrarily deep multilevel games) up to desired accuracy, and presented proofs for its feasibility and (under certain assumptions) optimality. We then used this algorithm to solve the multilevel problems we constructed and also solved a few problems from the literature for comparison, achieving satisfactory results in each case.

Future work in this direction would be to improve the convergence speed and accuracy of this algorithm. Furthermore, a desired generalization of this algorithm would be one that works with multiple leaders and multiple followers at all levels (or the so-called Multilevel Decentralized Problem). This would enable us to solve a wider variety of interesting problems that involve numerous stakeholders with varying levels of power amongst themselves. And finally, for applications where an exact solution is required, we want to explore methods to obtain them by leveraging the approximate solution provided by MCMO.

[Anandalingam, 1988]
authorAnandalingam, G., year1988.
titleA mathematical programming model of decentralized multi-level systems.
journalJournal of the Operational Research Society volume39, pages1021–1033.

[Bard, 1984]
authorBard, J.F., year1984.
titleAn investigation of the linear three level programming problem.
journalIEEE Transactions on Systems, Man, and Cybernetics volumeSMC-14, pages711–717.

[Bard and Moore, 1990]
authorBard, J.F., authorMoore, J.T., year1990.
titleA branch and bound algorithm for the bilevel programming problem.
journalSIAM Journal on Scientific and Statistical Computing volume11, pages281–292.

[Ben-Ayed and Blair, 1990]
authorBen-Ayed, O., authorBlair, C.E., year1990.
titleComputational difficulties of bilevel linear programming.
journalOperations Research volume38, pages556–560.

[Benson, 1989]
authorBenson, H.P., year1989.
titleOn the structure and properties of a linear multilevel programming problem.
journalJournal of Optimization Theory and Applications volume60, pages353–373.

[Bialas and Karwan, 1984]
authorBialas, W.F., authorKarwan, M.H., year1984.
titleTwo-level linear programming.
journalManagement science volume30, pages1004–1020.

[Blair, 1992]
authorBlair, C., year1992.
titleThe computational complexity of multi-level linear programs.
journalAnnals of Operations Research volume34.

[Cassidy et al., 1971]
authorCassidy, R.G., authorKirby, M.J.L., authorRaike, W.M., year1971.
titleEfficient distribution of resources through three levels of government.
journalManagement Science volume17, pages462–473.

[Dempe, 2020]
authorDempe, S., year2020.
titleBilevel optimization: Theory, algorithms, applications and a bibliography.

[Fard and Hajiaghaei-Keshteli, 2018]
authorFard, A.M.F., authorHajiaghaei-Keshteli, M., year2018.
titleA tri-level location-allocation model for forward/reverse supply chain.
journalAppl. Soft Comput. volume62, pages328–346.

[Fard et al., 2018]
authorFard, A.M.F., authorHajiaghaei-Keshteli, M., authorMirjalili, S.M., year2018.
titleHybrid optimizers to solve a tri-level programming model for a tire closed-loop supply chain network design problem.
journalAppl. Soft Comput. volume70, pages701–722.

[Han et al., 2015]
authorHan, J., authorLu, J., authorHu, Y., authorZhang, G., year2015.
titleTri-level decision-making with multiple followers: Model, algorithm and case study.
journalInformation Sciences volume311, pages182–204.

[Han et al., 2017]
authorHan, J., authorLu, J., authorZhang, G., year2017.
titleTri-level decision-making for decentralized vendor-managed inventory.
journalInf. Sci. volume421, pages85–103.

[Han et al., 2016]
authorHan, J., authorZhang, G., authorHu, Y., authorLu, J., year2016.
titleA solution to bi/tri-level programming problems using particle swarm optimization.
journalInf. Sci. volume370-371, pages519–537.

[Júdice and Faustino, 1992]
authorJúdice, J.J., authorFaustino, A.M., year1992.
titleA sequential lcp method for bilevel linear programming.
journalAnnals of Operations Research volume34, pages89–106.

[Labbé et al., 1998]
authorLabbé, M., authorMarcotte, P., authorSavard, G., year1998.
titleA bilevel model of taxation and its application to optimal highway pricing.
journalManagement science volume44, pages1608–1622.

[Lai, 1996]
authorLai, Y.J., year1996.
titleHierarchical optimization: A satisfactory solution.
journalFuzzy Sets Syst. volume77, pages321–335.

[Laine et al., 2023]
authorLaine, F., authorFridovich-Keil, D., authorChiu, C.Y., authorTomlin, C., year2023.
titleThe computation of approximate generalized feedback nash equilibria.
journalSIAM Journal on Optimization volume33, pages294–318.

[Liu, 1998]
authorLiu, B., year1998.
titleStackelberg-nash equilibrium for multilevel programming with multiple followers using genetic algorithms.
journalComputers & Mathematics with Applications volume36, pages79–89.

[Lu et al., 2016]
authorLu, J., authorHan, J., authorHu, Y., authorZhang, G., year2016.
titleMultilevel decision-making: A survey.
journalInf. Sci. volume346-347, pages463–487.

[Luo et al., 2020]
authorLuo, X., authorLiu, Y., authorLiu, J., authorLiu, X., year2020.
titleEnergy scheduling for a three-level integrated energy system based on energy hub models: A hierarchical stackelberg game approach.
journalSustainable Cities and Society volume52, pages101814.

[Oduguwa and Roy, 2002]
authorOduguwa, V., authorRoy, R., year2002.
titleBi-level optimisation using genetic algorithm, in: booktitleProceedings 2002 IEEE International Conference on Artificial Intelligence Systems (ICAIS 2002), organizationIEEE. pp. pages322–327.

[Pramanik and Roy, 2007]
authorPramanik, S., authorRoy, T.K., year2007.
titleFuzzy goal programming approach to multilevel programming problems.
journalEur. J. Oper. Res. volume176, pages1151–1166.

[Shih et al., 1996]
authorShih, H.S., authorLai, Y.J., authorLee, E.S., year1996.
titleFuzzy approach for multi-level programming problems.
journalComput. Oper. Res. volume23, pages73–91.

[Sinha, 2003]
authorSinha, S., year2003.
titleFuzzy programming approach to multi-level programming problems.
journalFuzzy sets and systems volume136, pages189–202.

[Tian et al., 2019]
authorTian, M., authorCui, M., authorDong, Z., authorWang, X., authorYin, S., authorZhao, L., year2019.
titleMultilevel programming-based coordinated cyber physical attacks and countermeasures in smart grid.
journalIEEE Access volume7, pages9836–9847.

[Tilahun, 2019]
authorTilahun, S.L., year2019.
titleFeasibility reduction approach for hierarchical decision making with multiple objectives.
journalOperations Research Perspectives .

[Tilahun et al., 2012]
authorTilahun, S.L., authorKassa, S.M., authorOng, H.C., year2012.
titleA new algorithm for multilevel optimization problems using evolutionary strategy, inspired by natural adaptation, in: booktitlePRICAI 2012: Trends in Artificial Intelligence: 12th Pacific Rim International Conference on Artificial Intelligence, Kuching, Malaysia, September 3-7, 2012. Proceedings 12, organizationSpringer. pp. pages577–588.

[Ue-Pyng and Bialas, 1986]
authorUe-Pyng, W., authorBialas, W.F., year1986.
titleThe hybrid algorithm for solving the three-level linear programming problem.
journalComputers & operations research volume13, pages367–377.

[Wächter and Biegler, 2006]
authorWächter, A., authorBiegler, L.T., year2006.
titleOn the implementation of an interior-point filter line-search algorithm for large-scale nonlinear programming.
journalMathematical programming volume106, pages25–57.

[White, 1997]
authorWhite, D.J., year1997.
titlePenalty function approach to linear trilevel programming.
journalJournal of Optimization Theory and Applications volume93, pages183–197.

[Woldemariam and Kassa, 2015]
authorWoldemariam, A.T., authorKassa, S.M., year2015.
titleSystematic evolutionary algorithm for general multilevel stackelberg problems with bounded decision variables (seamsp).
journalAnnals of Operations Research volume229, pages771–790.

[Yao et al., 2007]
authorYao, Y., authorEdmunds, T., authorPapageorgiou, D., authorAlvarez, R., year2007.
titleTrilevel optimization in power network defense.
journalIEEE Transactions on Systems, Man, and Cybernetics, Part C (Applications and Reviews) volume37, pages712–718.

[Zhang et al., 2010]
authorZhang, G., authorLu, J., authorMontero, J., authorZeng, Y., year2010.
titleModel, solution concept, and kth-best algorithm for linear trilevel programming.
journalInf. Sci. volume180, pages481–492.

§ NESTED TOLL-SETTING PROBLEM

The third level can be reformulated as:

\begin{align}
\min_{p_2} ~&p_2 \cdot (p_2 + t_2) + \\&(1-p_1-p_2) \cdot (1-p_1-p_2+D)\\
&0\le p_2 \le 1-p_1
\end{align}

The unconstrained stationary point for this problem is when:

\begin{align}
\nonumber p_2 + p_2 + &t_2 - (1-p_1-p_2)\\&- (1-p_1-p_2+D) = 0\\
&4p_2 + 2p_1 +t_2-D-2=0\\
&p_2 = \frac{D+2-2p_1-t_2}{4}
\end{align}

From the theory of constrained minimization, the response of the third level can then be written as:

\begin{align}
\label{eq:p2response}
p_2(t_1, t_2) = \begin{cases}
0 & \text{if} ~ D+2-2p_1(t_1)\le t_2 \\
1-p_1(t_1) & \text{if} ~D-2+2p_1(t_1)\ge t_2 \\
\frac{2+D-2p_1(t_1)-t_2}{4} & \text{otherwise}
\end{cases}
\end{align}

Similarly, we can obtain the response of the second level as:

\begin{align}
\label{eq:p1response}
p_1(t_1) = \begin{cases}
0 & \text{if} ~ t_1 \ge 2 \\
1 & \text{if} ~ t_1 \le -2\\
\frac{1}{4}(2-t_1) & \text{otherwise}
\end{cases}
\end{align}

From equations <ref> and <ref>, we can define the following parameterized constraint sets:

\begin{align*}
\mathcal{C}_1(D) := \{ & p_1=0 \wedge t_1 \ge 2, \\
& p_1=1 \wedge t_1\le-2, \\
&p_1=\frac{1}{4}(2-t_1) \wedge -2< t_1< 2\}
\end{align*}

\begin{align*}
\mathcal{C}_2(D) := \{ &p_2=0 \wedge ~ D+2-2p_1\le t_2 \\
&p_2=1-p_1 \wedge ~D-2+2p_1(t_1)\ge t_2 \\
&p_2=\frac{2+D-2p_1(t_1)-t_2}{4} \wedge \\
&~~~~~~~D-2+2p_1(t_1)< t_2 < D+2-2p_1 \}
\end{align*}

The solution for the Nested toll-setting problem would then simply be:

\begin{align}
\label{finaleq}
\max_{t_1, t_2, p_1, p_2} &p_1 \cdot t_1 + p_2 \cdot t_2 \\
\nonumber &t_1 \ge 0, t_2 \ge 0\\
\nonumber &(t_1,t_2,p_1,p_2) \in \bigcup \mathcal{C}_1(D) \cdot \mathcal{C}_2(D)
\end{align}

Equation <ref> is a standard quadratic programming problem defined over a union of polyhedral regions. It can be solved for each of the polyhedral regions using a standard solver to obtain the optimum value for the problem as follows:

* For $D = 6$, the obtained maximum is $4$ for $p_1=0, t_1=14, p_2=1, t_2=4$. The results imply that the toll-setter benefits when no traffic takes the tolled road at station $T_1$, i.e., ($p_1=0$) and all traffic takes the tolled road at station $T_2$.
It can be seen from equation <ref>a that the same objective can be realized for any $t_1\ge2$ (as this makes $p_1=0$).

* For $D = -1.5$, obtained maximum is $0.25$ for $p_1=0.25, t_1=1, p_2=0, t_2=4.53$. For this case, the toll-setter has to obtain all income from station $T_1$ as no traffic will take station $T_2$ due to the incentive on the non-tolled road, i.e., $p_2=0$. Like before, from equation <ref>a, for the given values of $p_1$ and $D$, any $t_2\ge0$ is a solution, as this yields $p_2=0$ for the same objective.