A simplex-type algorithm for continuous linear programs with constant coefficients

TL;DR

This paper introduces a simplex-type algorithm for solving continuous linear programs with constant coefficients over a finite time horizon, extending classical LP concepts to measure spaces and ensuring finite-step optimal solutions.

Contribution

It develops a finite-step simplex-like algorithm for continuous linear programs in measure spaces, generalizing previous models and their duals.

Findings

Algorithm solves problems in finite steps

Extends LP theory to measure and function spaces

Ensures existence of optimal strongly dual solutions

Abstract

We consider continuous linear programs over a continuous finite time horizon $T$, with a constant coefficient matrix, linear right hand side functions and linear cost coefficient functions, where we search for optimal solutions in the space of measures or of functions of bounded variation. These models generalize the separated continuous linear programming models and their various duals, as formulated in the past by Anderson, by Pullan, and by Weiss. In previous papers we have shown that these problems possess optimal strongly dual solutions. We also have presented a detailed description of optimal solutions and have defined a combinatorial analogue to basic solutions of standard LP. In this paper we present an algorithm which solves this class of problems in a finite bounded number of steps, using an analogue of the simplex method, in the space of measures.

Figures

Equations

$$\max$$

$$U(0-)=0$$

$$\min$$

$$P(0-)=0$$

$$\max$$

$${\int }_0^{t}Gu(s)ds+Fx(t)\le \alpha +at,$$

$$\max$$

$$\displaystyle\left[\begin{array}[]{cc}2&1\\ 1&1\end{array}\right]\left[\begin{array}[]{c}U_{1}(t)\\ U_{2}(t)\end{array}\right]+\left[\begin{array}[]{c}x_{1}(t)\\ x_{2}(t)\end{array}\right]=\left[\begin{array}[]{c}4\\ 1\end{array}\right]+\left[\begin{array}[]{c}4\\ 3\end{array}\right]t,$$

$$\begin{array}[]{clll}\gamma=[-4;-2],\\ c=[3;2],\\ A=\left[\begin{array}[]{cc}2&1\\ 1&1\end{array}\right],&b=\left[\begin{array}[]{c}4\\ 3\end{array}\right],&\beta=\left[\begin{array}[]{c}4\\ 1\end{array}\right],&T=3.\end{array}$$

$$\begin{array}[]{ccl}{\cal{L}}(0)=&(\beta_{0},\gamma_{0},\lambda_{0},\mu_{0},T_{0})&=\left(\left[\begin{array}[]{c}2\\ 1\end{array}\right],\left[\begin{array}[]{c}-5\\ -5\end{array}\right],\left[\begin{array}[]{c}-1\\ -1\end{array}\right],\left[\begin{array}[]{c}1\\ 1\end{array}\right],1\right),\\ {\cal{L}}(1)=&(\beta,\gamma,\lambda,\mu,T)&=\left(\left[\begin{array}[]{c}4\\ 1\end{array}\right],\left[\begin{array}[]{c}1\\ 2\end{array}\right],\left[\begin{array}[]{c}0\\ 0\end{array}\right],\left[\begin{array}[]{c}0\\ 0\end{array}\right],3\right),\end{array}$$

$$\max$$

$$AU(t)+x(t)=\beta +bt,0\le t<T,$$

$$\min$$

$$A^{\text{mbox}\mathsf{T}}P(t)-q(t)=\gamma +ct,0\le t<T,$$

$${\int }_{0-}^{T}x(T-t{)}^{\text{mbox}\mathsf{T}}dP(t)={\int }_{0-}^{T}q(T-t{)}^{\text{mbox}\mathsf{T}}dU(t)=0,$$

$$\max$$

$$A\mathbf{u}+AU\le \beta +bT+\lambda ,$$

$$\begin{array}[]{rcl}&\max&c^{{\mbox{\tiny\bf\sf T}}}u\\ &\mbox{s.t.}&Au+\dot{x}=b\\ \mbox{Rates-LP$({\cal{K}}_{n},{\cal{J}}_{n})$}&&u_{j}\in\mathbb{Z}\text{ for }j\in{\cal{J}}_{n},\;u_{j}\in\mathbb{P}\text{ for }j\notin{\cal{J}}_{n}\\ &&\dot{x}_{k}\in\mathbb{U}\text{ for }k\in{\cal{K}}_{n},\;\dot{x}_{k}\in\mathbb{P}\text{ for }k\notin{\cal{K}}_{n}\end{array}$$

$$\begin{array}[]{rcl}&\min&b^{{\mbox{\tiny\bf\sf T}}}p\\ &\mbox{s.t.}&A^{{\mbox{\tiny\bf\sf T}}}p-\dot{q}=c\\ \mbox{Rates-LP${}^{*}({\cal{K}}_{n},{\cal{J}}_{n})$}&&p_{k}\in\mathbb{Z}\text{ for }k\in{\cal{K}}_{n},\;p_{k}\in\mathbb{P}\text{ for }k\notin{\cal{K}}_{n}\\ &&\dot{q}_{j}\in\mathbb{U}\text{ for }j\in{\cal{J}}_{n},\;\dot{q}_{j}\in\mathbb{P}\text{ for }j\notin{\cal{J}}_{n}.\end{array}$$

$$x_k(t_n)=x_k^0+\sum _{m=1}^n{\dot{x}}_k^m{\tau }_m=0,$$

$$\sum _{n=1}^N{\tau }_n=T.$$

$$\begin{array}[]{lllll}\mathbf{u}^{0}_{j}=0,&j\in{\cal{J}}_{0}&\qquad\qquad\mathbf{p}^{N}_{k}=0,&k\in{\cal{K}}_{N+1},\\ x^{0}_{k}=0,&k\not\in{\cal{K}}_{0},&\qquad\qquad q^{N}_{j}=0,&j\not\in{\cal{J}}_{N+1},\\ \mathbf{p}^{0}_{k}=0,&k\in{\cal{K}}_{0},&\qquad\qquad\mathbf{u}^{N}_{j}=0,&j\in{\cal{J}}_{N+1},\\ {\mathbf{q}}^{0}_{j}=0,&j\not\in{\cal{J}}_{0},&\qquad\qquad{\mathbf{x}}^{N}_{k}=0,&k\not\in{\cal{K}}_{N+1}.\end{array}$$

$$\begin{array}[]{llll}A\mathbf{u}^{0}&+x^{0}=\beta,&\quad A^{{\mbox{\tiny\bf\sf T}}}\mathbf{p}^{N}&-q^{N}=\gamma.\end{array}$$

$$A{\mathbf{u}}^N+{\mathbf{x}}^N-x^N=\lambda ,A^{\text{mbox}\mathsf{T}}{\mathbf{p}}^0-{\mathbf{q}}^0+q^0=\mu .$$

$$\begin{array}[]{ll}\displaystyle x_{k}^{n}=x^{0}_{k}+\sum_{m=1}^{n}\dot{x}^{m}_{k}\tau_{m},&k=1,\ldots,K,\;n=1,\ldots,N,\\ \displaystyle q_{j}^{n}=q^{N}_{j}+\sum_{m={n+1}}^{N}\dot{q}^{m}_{j}\tau_{m},&j=1,\ldots,J,\;n=1,\ldots,N.\end{array}$$

$$MH=\mathbf{R}.$$

$$\begin{array}[]{ccccccccccccc}H=[&\mathbf{u}^{0}&x^{0}&\mathbf{p}^{N}&q^{N}&\tau&x^{1}\dots x^{N}&q^{0}\dots q^{N-1}&\mathbf{u}^{N}&{\mathbf{x}}^{N}&\mathbf{p}^{0}&{\mathbf{q}}^{0}&]^{{\mbox{\tiny\bf\sf T}}}\end{array}$$

$$\rho\to{\mathbf{R}}:{\mathbf{R}}=\left[\begin{array}[]{cccccccccccccccc}\beta^{{\mbox{\tiny\bf\sf T}}}&0&0&\gamma^{{\mbox{\tiny\bf\sf T}}}&0&0&0&T&0&0&\lambda^{{\mbox{\tiny\bf\sf T}}}&0&0&\mu^{{\mbox{\tiny\bf\sf T}}}&0&0\end{array}\right]^{{\mbox{\tiny\bf\sf T}}}$$

$$x_k^{\bullet }=\underset{0\le t<T}{\min }{x}_k(t),{\tilde{x}}_k=x_k^0-x_k^{\bullet },k=1,\dots ,K,$$

$$q_j^{\bullet }=\underset{0<t\le T}{\min }{q}_j(t),{\tilde{q}}_j=q_j^N-q_j^{\bullet },j=1,\dots ,J.$$

$$\tilde{U}={\int }_0^{T}u(t)dt,\tilde{P}={\int }_0^{T}p(t)dt,$$

$$\max$$

$$(x^0,q^N)$$

$$\min$$

$${}^{\ast }(x^0,q^N)$$

Full text

1 Introduction

This paper presents a finite, exact, simplex-type algorithm for the solution of the symmetric pair of dual continuous linear programs of the form:

[TABLE]

[TABLE]

Here $A$ is a $K\times J$ constant matrix, $\beta,b,\gamma,c$ are constant vectors of corresponding dimensions, the integrals are Lebesgue-Stieltjes, and include the jumps at 0 in $U,P$. The unknowns are vectors of cumulative control functions $U,P$ and vectors of non-negative slack or state functions $x,q$, over the time horizon $[0,T]$. It is convenient to think of dual time as running backwards, so that $P(T-t)$ is the vector of dual variables that correspond to the constraints of (1) at time $t$, and $U(t)$ correspond to the constraints of (1) at time $T-t$. The special feature here is that the controls are measures, or equivalently functions of bounded variation. We refer to the coefficients of the matrix $A$ as the structure parameters, to $b,\,c$ as the rate parameters, and to $\beta,\gamma,T$ as the boundary parameters of the M-CLP problem.

This paper continues research described in [36, 30, 31]. The algorithm for the solution of M-CLP builds on and extends the algorithm [36] for solution of SCLP, and is somewhat similar in principle. While some SCLP problems cannot be solved with the algorithm of [36], every SCLP problem can be solved as a special case of an M-CLP problem, by the algorithm presented in this paper. The problem (1) is solved parametrically. It starts from an artificial set of boundary parameters $T_{0},\beta_{0},\gamma_{0}$ that have a simple optimal solution, and moves in a finite bounded number of steps along a parametric straight line ${\cal{L}}(\theta),\,0\leq\theta\leq 1$, to the boundary parameters $\beta,\gamma,T$ of the original problem, where in each step an M-CLP pivot is performed.

At this stage it is too early to estimate the efficiency of this algorithm. Like standard LP simplex its worst case performance is exponential. It is however worth mentioning here, that we are currently in the process of re-programming the SCLP algorithm, beyond the pilot inefficient MATLAB implementation that was reported in [36]. While we have not yet stabilized the improved code, we made some improvements to the MATLAB current code. In comparing our current implementation of the simplex based algorithm [36], against an optimized uniform discretization solved by CPLEX we have found that the simplex based algorithm can outperform CPLEX both in a run-time (more then 10 times faster for some problems) and in the quality of the obtained solutions. Moreover we have found that the simplex based algorithm is able to solve very large problems (10000 controls and 10000 constraints), and these problems can only be solved by CPLEX with very coarse discretization. Based on these encouraging results, we expect that the M-CLP simplex based algorithm described in this paper, will be competitive with the discretization approach.

1.1 Background and motivation

These problems belong to the area of continuous linear programming (CLP), as pioneered by Bellman in 1953, and further discussed and investigated in [14, 20, 32, 33]. Bellman, was motivated by economic problems, in particular the optimization of continuous time Leontief input output economics, which also motivated Dantzig [9]. No direct solution method for these continuous time problems has emerged, but it was immediately observed that by discretizing time these problems could be approximated by standard linear programs. Among others, discrete time Leontief models motivated the Dantzig-Wolfe decomposition method [10]. In 1978 Anderson [1] formulated the less general separated continuous linear programming problem (SCLP), motivated by some job shop scheduling problems. Further work by Anderson and co-workers included [3] and a book [2], and much progress was achieved by Pullan [23, 24, 25, 26, 27], who formulated a non-symmetric dual, proved strong duality, and devised an algorithm based on a sequence of discretizations that converged to optimality. SCLP found several important applications, including beside job-shop scheduling also control of queueing networks [21], quickest transshipment problem[15], and general network flows over time [13, 12]. However, the solution methods for these problems remained by discretization of time. Several ingenious methods for efficient approximations to solutions of SCLP have appeared: Fleischer-Sethuraman [11] use geometrically scaled discretization, Bertsimas-Luo [6] use a non-convex quadratic programming formulation, and Bampou-Kuhn [4] use polynomial approximations.

The first algorithm to provide an exact finite solution method for SCLP was devised by Weiss [36] who considered the following version of SCLP:

[TABLE]

where $G,\,H,\,F$ are fixed matrices, $\alpha,\,a,\,\gamma,\,c,\,b,\,\,d$ are fixed vectors, and the unknowns are bounded measurable controls $u_{j}(t)$ and absolutely continuous states $x_{k}(t)$. For this problem Weiss has formulated a symmetric dual problem with bounded measurable controls $p_{k}(t)$ and states $q_{j}(t)$, found a combinatorial description of the extreme points of the feasible regions of the primal and dual, and devised a simplex type algorithm that solves these problems in a finite bounded number of steps, all of this under some sufficient conditions. We note that many of the applications of SCLP can indeed be formulated as (1.1).

The results of [36] left a large gap in the theory, since there is no strong duality between the SCLP problem (1.1) and its symmetric dual. Therefore it is possible that SCLP (1.1), or its dual, or both, have no optimal solution even if both are feasible. In any of these cases, the algorithm devised by Weiss fails to work.

The main idea in [30, 31] and in the current paper, is to formulate the problems M-CLP/M-CLP∗ (1), (1) in the space of measures, or equivalently in the space of functions of bounded variation. Already Pullan’s dual was formulated in the space of measures, but it was not a symmetric dual, and could not be used in a simplex type algorithm. The only other paper that we know of that looks for solutions in the space of measures is Papageorgiou [22]. The formulation of M-CLP/M-CLP∗ as (1), (1), with the unknowns as functions of bounded variation, turns out to achieve complete analogy with standard, finite dimensional linear programs (LP). It is shown in [30, 31] that strong duality holds, with a complete combinatorial description of the finite number of extreme points of the feasible region. The last step in the analogy with LP is a finite simplex type algorithm derived in this paper. Here M-CLP stands for CLP in the space of measures.

1.2 Overview and contribution

At this point it is useful to consider the differences between SCLP (1.1), and M-CLP (1). Optimal solutions of SCLP have piecewise constant control rates $u_{j}(t)$ and continuous piecewise linear states $x_{k}(t)$, with complementary slack dual controls $p_{k}(T-t)$ and dual states $q_{j}(T-t)$. Hence the solution can be summarized by the breakpoints $0=t_{0}<\cdots<t_{N}=T$, by the values of $u_{j}(t)$ and of the derivatives $\dot{x}_{k}(t)$ and their duals in each interval, and by the boundary values of the primal $x_{k}(0)$, and dual $q_{j}(0)$. These are determined as follows: The boundary values $x_{k}(0)$ and $q_{j}(0)$ are determined by two standard LPs, with coefficient matrices $F$ and $H$ respectively. The rates $u_{j},\dot{x}_{k},p_{k},\dot{q}_{j}$ in each interval are a pair of complementary slack basic solution of an LP (referred to as the Rates-LP) with coefficient matrix \left[\begin{array}[]{cc}G&F\\ H&0\end{array}\right] (replaced by $A$ in M-CLP), and each of the breakpoints is determined by an equation of the form $x_{k}(t_{n})=0$ or of the form $q_{j}(T-t_{n})=0$. It follows from this structure that the optimal solution to an SCLP is determined by an optimal sequence of adjacent (differ by a single LP pivot) bases $B_{1},\ldots,B_{N}$ of the Rates-LP, and given this sequence the optimal solution can be computed directly. The SCLP algorithm constructs this optimal sequence of bases.

The algorithm solves SCLP parametrically, over the parametric line of time horizons moving from [math] to $T$. This is similar to the parametric self dual simplex algorithm, also known as Lemke’s algorithm [17, 34]. As a first step, the vectors of boundary values of the optimal solution, $x(0),\,q(0)$ are calculated, and these are valid and remain fixed for all time horizons. The algorithm moves from 0 to $T$ in a finite bounded number of steps, at each of these steps an SCLP pivot is performed, to update the sequence of bases. In this SCLP pivot some bases are deleted and some bases are added to the sequence of bases.

The optimal solution of M-CLP in the open interval $t\in(0,T)$ is likewise given by continuous piecewise linear $U_{j}(t)$ with piecewise constant derivatives $u_{j}(t)$ and continuous piecewise linear $x_{k}(t)$ with complementary slack duals $P_{k}(T-t)$ and $q_{j}(T-t)$, which are again determined by a sequence of adjacent bases $B_{1},\ldots,B_{N}$ of a Rates-LP with coefficient matrix $A$. In addition the solution includes impulse controls at 0 and $T$. Thus, the boundary values of optimal M-CLP/M-CLP∗ solutions consist of primal and dual impulse controls $U(0)$, $U(T)-U(T-)$, and $P(0)$, $P(T)-P(T-)$, and of the jumps in the states, $x_{k}(0)$, $x_{k}(T)-x_{k}(T-)$, $q_{j}(0)$, $q_{j}(T)-q_{j}(T-)$ (here we define by convention $x(0-)=0,\,q(0-)=0$). These results were derived in [30, 31]. The solution of M-CLP can be decomposed into a solution for an Internal-SCLP problem, solved by a sequence of adjacent bases $B_{1},\ldots,B_{N}$ of $[A\;I]$ (which can be solved by the algorithm of [36]), and a pair of Boundary-LP problems. One of the boundary LPs determines the optimal primal boundary values at [math] and $T$, the other determines the optimal dual boundary values at [math] and $T$. The two boundary problems are of dimension $2K\times 2J$ each. They share the same (transposed) matrix of coefficients, and their solutions are complementary slack but not dual.

The algorithm for M-CLP uses a more complicated parametric line ${\cal{L}}(\theta)$, $0\leq\theta\leq 1$, that moves from initial values $\beta_{0},\gamma_{0},\lambda_{0},\mu_{0},T_{0}$ to the final values $\beta,\gamma,0,0,T$ ($\lambda,\,\mu$ are additional vectors of boundary parameters needed here, but they disappear in the final solution, they will be explained in Section 3 and Appendix A.2). The solution again moves along the parametric line in a finite number of steps, where at each step an M-CLP pivot is performed. This M-CLP pivot may be an internal pivot on the Internal-SCLP or it may be a boundary pivot, that involves the Boundary-LPs.

Two issues had to be resolved in the derivation of the M-CLP algorithm. The decomposition into Boundary-LPs and Internal-SCLP is quite subtle, it will be derived in Section 4. Also, a significant difference between the SCLP and the M-CLP algorithm is that for SCLP the equations for the boundary values, for the time intervals, and for the values of the state variables at the breakpoints are decoupled, while in the M-CLP problem these equations are coupled. These equations need to be solved to obtain the next point of the parametric line at which a constraint becomes tight and the next iteration is performed, as explained in Section 8.

1.3 Applications and extensions

We now list some potential applications and possible extensions of M-CLP.

Solution of SCLP porblems

By formulating SCLP as M-CLP, problems that cannot be solved by the algorithm of [36] can now be solved.

Sensitivity analysis:

The fact that we obtain an exact solution, with a complete detailed structure of the solution, enables us to do sensitivity analysis similar to what is done with standard LP. For example, we can see how the solution may change if we change the time horizon, or if we move some of the time points $t_{n}$. This is much harder to do if the problem is solved by discretization.

Rolling horizon and model predictive control

In many applications of optimization over time, problems are solved for some given time horizon $T$, and are then re-solved repeatedly, for the same length of time horizon, after every $d=T/m$ time units, with updated data. For M-CLP that would mean updated $\beta,\gamma$. By using the simplex type algorithm for M-CLP, this can be done very efficiently, since updating requires only a fraction of the calculations required to solve the complete problem. This efficiency is lost if one uses discretization.

Continuous time Leontief economic model

We can formulate and solve Leontief economic model problems in continuous time, using an M-CLP formulation, if the target capacity function is given in advance (cf. [29] Section 1.3). If one wishes also to optimize the capacity this becomes a CLP of the type formulated by Bellman [5] which is not an M-CLP.

SCLP and M-CLP with piecewise constant data

The M-CLP algorithm can directly be generalized to solve problems with piecewise constant $b,c,\beta,\gamma$. Solutions of such problems will have impulse controls at points at which the data changes.

Fractional continuous linear programs

Problems with linear constraints and an objective function that is the ratio of two linear terms are known as fractional linear programs, and can be solved by a simplex algorithm. A continuous time version of this was studied by Wen [37, 38], who analyzed their properties and suggested an approximation algorithm that included bounding the objective and discretizing time. It is possible to formulate such problems as M-CLP, and use our algorithm to obtain the exact solution in a finite number of steps.

Continuous linear complementarity problems

Linear complementarity problems (LCP) have been studied since the 60’s (cf. comprehensive book of Cottle, Pang and Stone [8]). While general LCP are NP-hard, several types of LCP’s can be solved by a simplex type algorithm developed by Lemke and Howson [16, 17, 18]. These include convex quadratic programming problems, and bi-matrix games. A continuous time version of LCP can be formulated, and a simplex type algorithm that is similar to our M-CLP algorithm can be used to solve these problems.

More general right hand side and objective functions

Pullan [24] has shown that if the right hand side and objective functions are piecewise analytic, and the feasible region is bounded, SCLP has a piecewise analytic optimal solution. We hope to obtain similar results for M-CLP, and extend the algorithm to such problems. It is reasonable to conjecture that optimal solutions will again consist of a finite partition of the time horizon, where each interval is associated with a basis of $[A\;I]$, and the control rates and rates of the states are analytic functions in each interval. It seems that one can then extend our algorithm to such r.h.s. and objective, with the added difficulty that along the parametric line, constraints become tight not just at the breakpoints, but also at internal points of intervals, at times at which the derivative of a state variable is 0 (i.e. at local minima).

1.4 Outline or the rest of the paper

The rest of the paper is constructed as follows: In Section 2 we give an example of an M-CLP problem solved by our algorithm. It is the smallest possible example, but was chosen to have enough iterations to illustrate most features of the algorithm. This should provide some intuition, before we introduce the actual algorithm in Section 8. We return to this example in Section 9.

In the next part of the paper we extend some of the results of [31]: we give a succinct description of the solution (Section 3), then describe a decomposition of M-CLP and its dual to an Internal-SCLP problem and a pair of Boundary-LP problems (Section 4), and follow with a discussion of validity regions (Section 5).

To construct the algorithm, we discuss in Section 6 what happens on the boundaries of validity regions. When ${\cal{L}}(\theta)$ approaches the boundary of a validity region, at the boundary, some of the values of the impulse controls or the state variables or the interval lengths, which are positive in the interior of the validity region, reach 0. As a result, the set of positive boundary values and sequence of internal bases is no longer optimal outside the boundary. We call this a collision, and we classify all possible types of collisions.

In Section 7 we discuss M-CLP pivots. We consider two neighboring validity regions ${\cal{V}},{\cal{W}}$, and given the type of collision as ${\cal{L}}(\theta)$ approaches the boundary of ${\cal{V}}$, we show how to identify the collision as ${\cal{L}}(\theta)$ approaches the point of collision from the other direction, i.e. from within ${\cal{W}}$. This enables us to identify the set of positive boundary values and the sequence of internal bases that are optimal in ${\cal{W}}$, which completes the pivot.

Based on the pivot operation we then describe the algorithm in Section 8, where we show how to determine the collision points and perform the sequence of pivots. In Section 9 we return to the example of Section 2, and re-interpret it in terms of the algorithm.

The algorithm of Section 8 requires us to make some strong simplifying assumptions in order for it to function. However, in all cases these assumptions can be enforced by perturbation of the parametric line. We outline the necessary perturbation procedures in Appendix A.1.

For ease of reading and continuity we moved some of the proofs of Sections 5, 7 and 8 to Appendix A.3, A.4, A.5 respectively.

Notes on presentation: Throughout the paper we quote results from [30, 31]. We will refer to results from the duality paper [30] by appending D and from the structure paper [31] by appending S (e.g. D5.3, or (S4.2)). There is complete symmetry between M-CLP and M-CLP∗. Therefore when we formulate results for one or both of them, we will give a proof only for one of them.

2 An Illustrative Example

We now apply the algorithm to solve the following problem:

[TABLE]

As can be seen by comparing with (1), the data of this problem is:

[TABLE]

Figure 1 gives a graphical representation of the complete solution of this problem, as well as of its dual. The challenge in presenting this solution in a graph is that we wish to show 8 different function variables in one picture, in addition to impulse controls at 0 and $T$. We also wish to distinguish the different types of variables and illustrate the duality. To do so we use the horizontal axis of length $T$ to present primal time from left to right and dual time in the reversed direction from right to left. We present primal states with the complementary slack dual control rates in the upper quadrant, and the dual states with the complementary slack primal control rates in the lower quadrant. To be precise, we plot all the the non-zero values of: $x_{k}(t)$, $p_{k}(T-t)$, $-q_{j}(T-t)$, and $-u_{j}(t)$. We do not plot 0 values, and because basic solutions of the $2\times 2$ Rates-LP have only 2 non-zero primal values and 2 non-zero dual values, we only need to plot two positive and two negative values for each $0<t<T$. It is easy to distinguish between state variables ($x_{k},\,q_{j}$) which are continuous piecewise linear, and control rates (derivatives of the cumulative controls, $u_{j},\,p_{k}$) which are piecewise constant.

The impulse controls are shown as vertical arrows, next to the origin and next to the horizontal coordinate $T$. Next to the origin we plot $-U_{j}(0)$ and $P_{k}(T)-P_{k}(T-)$. Next to $T$ we plot $P_{k}(0)$ and $-(U_{j}(T)-U_{j}(T-)$. Circles are used at the origin for jumps in $q_{j}(T)$, and at $T$ for jumps in $x_{k}(T)$.

We use different colors to refer to the primal and dual states, and for each state we use the same color also for the complementary control. We use the same colors in all the following figures.

In this figure we can see that the solution consists of 3 intervals, and a single impulse control at 0 for the primal and another at 0 for the dual (in the reversed time). Circles denote the values of the state variables at $T$, which might be different from the limiting values at $T{-}$.

To apply the algorithm we will solve the problem along a parametric line, that will move along the line of boundary parameters ${\cal{L}}(\theta)$ from $\theta=0$ to $\theta=1$, where

[TABLE]

The artificial boundary parameters $\lambda,\mu$ are used to avoid degeneracies along ${\cal{L}}(\theta)$, and disappear at $\theta=1$. They will be explained in Appendix A.2.

The algorithm will solve this problem in 7 iterations, the problem was deliberately chosen to have many iterations, so as to demonstrate various features of the algorithm. The iterations happen at $\theta$ values $0,\frac{1}{11},\frac{1}{6},\frac{2}{9},\frac{4}{9},\frac{5}{9},\frac{13}{19},1$. In each of the intervals between these $\theta$s, the optimal solution maintains the same set of positive boundary values, the same number of breakpoints in the solution, and the same rates of controls and slopes of the states between breakpoints.

As $\theta$ changes, the solution changes continuously in an affine way, until a boundary constraint or an internal constraint becomes tight. At this $\theta$ the edge of the validity region is reached and the pivot of the iteration happens. If a boundary value shrinks to zero then we perform a boundary pivot. Otherwise, if some intervals shrink to zero, or some primal or dual state variable, at a breakpoint, shrinks to zero, we perform an Internal-SCLP pivot. We call this an M-CLP pivot.

The initial values ${\cal{L}}(0)$ are chosen so that the solution has a single interval, with 0 controls and states equal to the right hand sides: $x(t)=\beta(\theta)+bt$, q(t)=-\big{(}\gamma(\theta)+ct\big{)}, $0<t<T(\theta)$, for small $\theta$.

In the following paragraphs we present illustrations of the optimal solution at each iteration, showing the solution at the boundary and the middle of the interval of $\theta$s, and we briefly explain what happens at the pivot. We return to this example in Section 9, with additional details, after the algorithm is derived.

Iteration 1: Starting at $\theta=0$ the solution is valid until $\theta=1/11$, when $q_{1}(T)$ shrinks to 0 (at the left side of the figure, recall dual time runs right to left).

A boundary pivot is performed and as a result a primal impulse control $U_{1}(0)$ at time 0 is introduced, and the boundary state $x_{2}(0)$ jumps down to 0. The boundary pivot also introduces a dual impulse control $P(T)-P(T-)$ that is still 0 but will increase in the next validity region.

Iteration 2: Starting at $\theta=1/11$ the solution is valid until $\theta=1/6$, when $q_{1}(T-)$ shrinks to zero. In the figure we show the solution also at $\theta=1/8$. Note that the impulse control $P(T)-P(T-)$ at the left side of the figure is now positive and growing.

An Internal-SCLP pivot is performed, in which a new interval is inserted at primal time 0 (dual time $T$), where the state variable $q_{1}$ shrunk to 0. Currently the new interval has length 0. It will increase with $\theta$.

Iteration 3: Starting at $\theta=1/6$ the solution now consists of two intervals, $0<t_{1}<T$. At $\theta=2/9$ the value of $x_{1}(t_{1})$ shrinks to 0. This happens at $t_{1}=2/9$, when the local minimum of $x_{1}(t)$ descends to 0. We also show the solution at $\theta=7/36$, where the local minimum of $x_{1}(t)$ is still $>0$. Note that the new interval has increased in length, notice also that there are positive primal and dual control rates in this interval.

An Internal-SCLP pivot is performed, and a new interval is inserted at $t_{1}$. Currently the new interval has length 0.

Iteration 4: Starting at $\theta=2/9$ the solution now consists of three intervals, $0<t_{1}<t_{2}<T$. At $\theta=4/9$ the value of $q_{2}(T-t_{1})$ shrinks to 0. We also show $\theta=3/9$.

An Internal-SCLP pivot is performed, and another new interval is inserted at $t_{1}$. It is currently of length 0.

Iteration 5: Starting at $\theta=4/9$ the solution now consists of four intervals, $0<t_{1}<t_{2}<t_{3}<T$. At $\theta=5/9$ the interval between $t_{2}$ and $t_{3}$ shrinks to 0. We also show $\theta=1/2$. Note that in the interval $(t_{1},t_{2})$ all the states are zero, and all four control rates are positive.

An Internal-SCLP pivot is performed, and the positive variables in the 3rd interval are changed. What happened was that for $4/9<\theta<5/9$, at $t_{2}$ we had $x_{2}(t_{2})=q_{2}(T-t_{2})=0$, and at $t_{3}$, $x_{1}(t_{3})=q_{1}(T-t_{3})=0$. As $\theta\to 5/9$, $t_{2}$ and $t_{3}$ moved toward each other until they met. After the pivot, with $\theta>5/9$ they continued to move in the same directions, so now the roles of $t_{2},t_{3}$ have switched, and we will have for $\theta>5/9$ that $x_{2}(t_{3})=q_{2}(T-t_{3})=0$ and $x_{1}(t_{2})=q_{1}(T-t_{2})=0$. So in this pivot, the interval $(t_{2},t_{3})$ shrunk to zero, and after the pivot it will increase again, but with different positive state variables.

Iteration 6: Starting at $\theta=5/9$ the solution now consists of four intervals, $0<t_{1}<t_{2}<t_{3}<T$, where the interval $(t_{2},t_{3})$ now has a different basic solution of the Rates-LP. At $\theta=13/19$ the boundary value of $x_{2}(T)$ shrinks to 0.

A boundary pivot is performed. A new impulse control $P_{2}(0)$ is introduced at dual time 0 (primal time $T$, i.e. at the right end of the illustration). This impulse control is currently still 0.

Iteration 7: Starting at $\theta=13/19$ the solution still consists of four intervals. At $\theta=1$, the interval between $t_{3}$ and $t_{4}$ shrinks to 0, and only three intervals remain. Also the impulse control $P_{2}(T)-P(T-)$, at the left end of the illustration, shrinks to 0.

The result at $\theta=1$ is the optimal solution to our problem, as shown previously in Figure 1.

3 Structure of the Solution

All the results in this section follow directly from [30, 31], and Appendix A.2. We therefore state them here briefly, and the reader is advised that the notation and terms presented here will be used for the rest of the paper.

We formulate a slightly more general M-CLP problem and its dual, which our algorithm will solve:

[TABLE]

[TABLE]

The new elements in this re-formulation are two additional vectors of boundary parameters: $\lambda\leq 0,\,\mu\geq 0$. The introduction of $\lambda,\,\mu$ here is done to avoid some inherent degeneracies that can inhibit our algorithm. We discuss the motivation and explain the sign restrictions on $\lambda,\mu$, and we prove that all the results of [30, 31] extend directly to the current formulation, in Appendix A.2.

We denote the complete set of boundary parameters by $\rho=(\beta,\gamma,\lambda,\mu,T)$. We will often make the following non-degeneracy assumption:

Assumption 1** (Non-Degeneracy).**

The vector $b$ is in general position to the matrix $\left[A\;I\right]$ (i.e. it is not a linear combination of any less than $K$ columns), and the vector $c$ is in general position to the matrix $\left[A^{{\mbox{\tiny\bf\sf T}}}\;I\right]$.

Optimal solutions of (3), (3) are strongly dual and satisfy the following complementary slackness condition:

[TABLE]

Based on ideas of Wang, Zhang and Yao [35], and similar to Theorem S3.1, a necessary and sufficient condition that (3) be feasible is that the following Test-LP is feasible

[TABLE]

An analogous condition applies for M-CLP∗ (3). Feasibility of (3), (3) implies existence of optimal complementary slack solutions with equal objective values.

Under the Non-Degeneracy Assumption 1 optimal solutions $U(t),P(t)$ have impulse controls $\mathbf{u}^{0}=U(0),\,\mathbf{u}^{N}=U(T)-U(T-),\,\mathbf{p}^{0}=P(T)-P(T-),\,\mathbf{p}^{N}=P(0)$ at [math] and $T$, piecewise constant control rates $u(t)=\frac{dU(t)}{dt},\,p(t)=\frac{dP(t)}{dt}$ at $0<t<T$, and continuous piecewise linear states $x(t)=\beta+bt-U(t),\,q(t)=\gamma+ct-P(t)$ at $0\leq t<T$, with possible discontinuities at $T$.

The following list adds to the description of optimal solutions.

The time horizon $[0,T]$ is partitioned by $0=t_{0}<t_{1}<\ldots<t_{N}=T$ which are the breakpoints in the rates $u,p$ and in the slopes of $x,q$.

• -

We denote the values of interval lengths by $\tau_{n}=t_{n}-t_{n-1},\,n=1,\ldots,N$.

• -

We denote the vectors of values of the states at the breakpoints by $x^{n}=x(t_{n}),\,n=0,\ldots,N-1$, and $q^{n}=q(T-t_{n})$, $n=1,\ldots,N$.

• -

Because there may be a discontinuity at $T$ we denote the values at $T$ itself by ${\mathbf{x}}^{N}=x(T)$, ${\mathbf{q}}^{0}=q(T)$, and let $x^{N}=x(T-)$, $q^{0}=q(T-)$ be the values of the limits as $t\nearrow T$.

• -

The constant slopes of the states and the constant values of the control rates for each interval are denoted $\dot{x}^{n}=\frac{dx(t)}{dt},\,u^{n}=u(t),\,t_{n-1}<t<t_{n}$ and $\dot{q}^{n}=\frac{dq(t)}{dt},\,p^{n}=p(t),\,T-t_{n}<t<T-t_{n-1}$.

The entire solution, $U(t),P(t),x(t),q(t)$ can be retrieved from the vectors of boundary values $\mathbf{u}^{0},\mathbf{u}^{N},\mathbf{p}^{0},\mathbf{p}^{N}$, $x^{0},{\mathbf{x}}^{N},{\mathbf{q}}^{0},q^{N}$, the vectors of rates $u^{n},\dot{x}^{n},p^{n},\dot{q}^{n},$ $n=1,\ldots,N$, and and the interval lengths $\tau_{n}\,n=1,\ldots,N$.

We now continue to describe properties of optimal solutions. The rates $u^{n},\dot{x}^{n},p^{n},\dot{q}^{n}$ are complementary slack basic solutions of the following dual pairs of Rates-LP problems

[TABLE]

[TABLE]

Here $\mathbb{P},\mathbb{Z},\mathbb{U}$ refer to sign restrictions on the variables: $\mathbb{P}$ means non-negative, $\mathbb{U}$ means unrestricted, $\mathbb{Z}$ means restricted to be 0. For $n=1,\ldots,N$, $B_{1},\ldots,B_{N}$ denote the optimal bases, and ${\cal{K}}_{n},{\cal{J}}_{n}$ are the indexes of the basic $\dot{x}_{k}^{n},\dot{q}_{j}^{n}$. By Assumption 1 all the bases are non-degenerate. We denote by $B^{*}_{n}$ the dual basis complementary to $B_{n}$.

Definition 3.1**.**

We say that a base sequence $\{{\cal{K}}_{n},{\cal{J}}_{n}\}_{n=0}^{N+1}$ is a proper base sequence if it satisfies:

• •

The bases $B_{1},\ldots,B_{N}$ are admissible in the sense that $u^{n},p^{n}\geq 0$.

• •

For each $n=1,\dots,N-1$ bases $B_{n},B_{n+1}$ are adjacent, in the sense that $B_{n}\to B_{n+1}$ involves a single pivot, with $v_{n}$ leaving the basis, and $w_{n}$ entering.

• •

The indexes of the non-zero boundary values of the primal and dual state variables at times 0 and $T$ denoted by $\{{\cal{K}}_{n},{\cal{J}}_{n}\}_{n\in\{0,N+1\}}$ satisfy a compatibility condition that ${\cal{K}}_{0}\subseteq{\cal{K}}_{1},\,{\cal{J}}_{N+1}\subseteq{\cal{J}}_{N}$, i.e. if $x_{k}(0)>0$ then $\dot{x}_{k}^{1}\neq 0$, and if $q_{j}(0)>0$ then $\dot{q}_{j}^{N}\neq 0$.

Given any proper base sequence, we can attempt to construct an optimal solution by solving the Rates-LP/LP∗ (24), (25) for the bases $B_{n},\,n=1,\ldots,N$ and then using these values to formulate and solve the following coupled linear equations, which determine the interval lengths and the boundary values.

The time interval equations:

[TABLE]

[TABLE]

The complementary slackness conditions:

[TABLE]

The first boundary equations:

[TABLE]

The second boundary equations:

[TABLE]

Where all the values of the state variables are obtained from:

[TABLE]

Theorem S3.1 in [31] states that for a proper base sequence, if the solution to (26)–(31) is non-negative then it is an optimal solution, and every feasible M-CLP/M-CLP∗ problem has an optimal solution of this form, i. e. one can construct a proper base sequence such that solution of (26)–(31) will be non-negative.

Given a proper base sequence, after solving (24), (25), we can formulate all the equations (26)–(31) as

[TABLE]

The matrix $M$ is a square matrix of dimension $(J+K)(N+4)+N$, defined in Section S3.5 111There are some typos in the definition of $M$ in [31]. They were corrected in [29].. The vector of unknowns (writing all vectors as row vectors) is:

[TABLE]

and the right hand side consists of the boundary parameters $\rho$ interspersed by properly dimensioned 0 vectors:

[TABLE]

The set of equations (32) are formulated in such a way that they force $2K+2J$ boundary variables, $N-1$ variables that are left-hand sides of (26) and $N(K+J)-\sum_{n=1}^{N}(|{\cal{K}}_{n}|+|{\cal{J}}_{n}|)$ of the $x_{k}^{n}\in\{n,k:\dot{x}^{n}_{k}=0\},\,q_{j}^{n}\in\{n,j:\dot{q}^{n+1}_{j}=0\}$ that correspond to non-basic variables in the Rates-LP, to be zero. We denote this set of zero-valued variables by $H_{\mathbb{Z}}$. The set of remaining, non-negative variables, is denoted by $H_{\mathbb{P}}$. We say that the solution $H$ is fully non-degenerate if all components of $H_{\mathbb{P}}$ are $>0$.

4 Decomposition

We introduce the following distinction between the state variables that are indexed by ${\cal{K}}_{0},{\cal{J}}_{N+1}$:

Definition 4.1**.**

Consider a fully non-degenerate solution, and the $x_{k}(t)$ with $k\in{\cal{K}}_{0}$. If for some $t_{n}<T$, $x_{k}(t_{n})=0$ then we say that $x_{k}$ is tied. On the other hand, if $x_{k}(t)>0$ for $0\leq t<T$ then we say that $x_{k}$ is free. We denote by ${\cal{K}}^{=}$ the indexes of the tied variables, and by ${\cal{K}}^{\uparrow\downarrow}$ the indexes of the free variables, so that ${\cal{K}}_{0}={\cal{K}}^{=}\cup{\cal{K}}^{\uparrow\downarrow}$. We define ${\cal{J}}^{=},\,{\cal{J}}^{\uparrow\downarrow}$ analogously for $q_{j}(t),j\in{\cal{J}}_{N+1}$.

The tied state variables appear both in the boundary equations (29) and in the time interval equations (31), and introduce coupling between these different sets of equations. The free state variables appear only in the boundary equations (29).

Given an optimal solution of an M-CLP/M-CLP∗ problem, with data $A,b,c,\beta,\gamma,\lambda,\mu,T$, we introduce some more notations for various quantities:

[TABLE]

Here $x_{k}^{\bullet}$ is the minimum of $x_{k}(t)$ over the range $0\leq t<T$ (excluding the value $x_{k}(T)={\mathbf{x}}^{N}$). The values of $\tilde{x}_{k}$ can also be expressed as ${\tilde{x}}_{k}=-\min_{0\leq n\leq N}\sum_{m=1}^{n}\dot{x}_{k}^{m}\tau_{m}$ (where empty sum equals 0). For $k\in{\cal{K}}^{=}$, $x_{k}^{\bullet}=0$ and $\tilde{x}_{k}=x_{k}^{0}$. For $k\in{\cal{K}}^{\uparrow\downarrow}$, $x_{k}^{\bullet}>0$, and $\tilde{x}_{k}=x_{k}^{0}-x_{k}^{\bullet}\geq 0$. For $k\not\in{\cal{K}}_{0}$, we have $x_{k}^{0}=x_{k}^{\bullet}=\tilde{x}_{k}=0$. Analogous statements hold for $q_{j}(t)$.

We also introduce the following notation for the vectors of total cumulative controls excluding the impulse controls:

[TABLE]

We refer to ${\tilde{U}},{\tilde{P}},{\tilde{x}},{\tilde{q}}$ as decomposition parameters.

Proposition 4.2**.**

Under Non-Degeneracy Assumption 1 M-CLP/M-CLP∗ possess a unique set of decomposition parameters ${\tilde{U}},{\tilde{P}},{\tilde{x}},{\tilde{q}}$. They are all of them affine functions of elements of $\rho$.

Proof.

Uniqueness follows directly from Theorem D5.5(iii), which states that under Assumption 1 $u(t),p(t)$ are uniquely determined, and as a result so are $\dot{x}(t),\dot{q}(t)$. Proof that they are affine functions of $\rho$ follows similar to the proof of Corollary S3.8.   ∎

We are now ready to define a decomposition of M-CLP/M-CLP∗ to two parametric families of problems. The first consists of SCLP/SCLP∗ problems and concerns the solution in the interior of the time horizon. The second consists of pairs of Boundary-LP problems, and concerns the solution on the boundaries of the time horizon.

Internal-SCLP Problems:

Consider M-CLP/M-CLP∗ problems with data $A,b,c$. We define a parametric family SCLP/

SCLP${}^{*}(x^{0},q^{N})$ with non-negative parameters $x^{0},q^{N}$

[TABLE]

[TABLE]

These problem are a special case of separated continuous linear program (SCLP) of the form discussed in [36]. We call these problems Internal-SCLP/ SCLP∗.

Boundary Problems

Consider an M-CLP/M-CLP∗ problem with data $A,b,c,\beta,\gamma,T,\lambda,\mu$. We define, similar to S15, S16 in [31], a family of pairs of LP problems, for parameter vectors $\tilde{u},\tilde{x},{\tilde{U}},{\tilde{P}}$:

[TABLE]

[TABLE]

Note, that Boundary-LP$({\tilde{U}},\tilde{x})$ and Boundary-LP${}^{*}(\tilde{P},\tilde{q})$ share the same matrix of constraint coefficients (transposed), but they are not dual to each other. Note also that they use disjoint sets of decomposition parameters (${\tilde{U}},{\tilde{x}}$ for Boundary-LP and $\tilde{P},\tilde{q}$ for Boundary-LP∗).

Theorem 4.3**.**

(i) An optimal solution of M-CLP/M-CLP∗ can be decomposed to optimal solutions of Internal-SCLP/SCLP${}^{*}(x^{0},q^{N})$ and feasible solutions of Boundary-LP$({\tilde{U}},{\tilde{x}})$ and Boundary-LP${}^{*}({\tilde{P}},{\tilde{q}})$ satisfying:

[TABLE]

(ii) Conversely, a combination of optimal solutions of Internal-SCLP/ SCLP${}^{*}(x^{0},q^{N})$, and of feasible solutions of Boundary-LP$({\tilde{U}},{\tilde{x}})$ and Boundary-LP${}^{*}({\tilde{P}},{\tilde{q}})$, that satisfy (41) can be composed into an optimal solution of M-CLP/M-CLP∗.

(iii) The solutions of Boundary-LP$({\tilde{U}},{\tilde{x}})$ and Boundary-LP${}^{*}({\tilde{P}},{\tilde{q}})$ considered in (i) and (ii) are optimal.

Proof.

(i) Let $u(t),x(t),p(t),q(t),\mathbf{u}^{0},\mathbf{u}^{N},x^{0},{\mathbf{x}}^{N},\mathbf{p}^{N},\mathbf{p}^{0},q^{N},{\mathbf{q}}^{0}$ be an optimal solution of M-CLP/M-CLP∗ as described in Theorem S3.1. Then it is immediate to see that the same $u(t),x(t),p(t),q(t)$ solve the Internal-SCLP/SCLP∗ $(x^{0},q^{N})$.

Furthermore, one can see that $\mathbf{u}^{0},\mathbf{u}^{N},x^{\bullet}=x^{0}-{\tilde{x}},{\mathbf{x}}^{N}$ is a feasible solution of Boundary-LP$({\tilde{U}},{\tilde{x}})$ and $\mathbf{p}^{N},\mathbf{p}^{0},q^{\bullet}=q^{N}-{\tilde{q}},{\mathbf{q}}^{0}$ is a feasible solution of Boundary-LP${}^{*}({\tilde{P}},{\tilde{q}})$.

(ii) On the other hand, optimal solution of SCLP/SCLP∗ provides $u(t),$ $p(t),x(t),q(t)$ and solutions of Boundary-LP/LP∗ provides $\mathbf{u}^{0},\mathbf{u}^{N},x^{\bullet},{\mathbf{x}}^{N},\mathbf{p}^{0},$ $\mathbf{p}^{N},q^{\bullet},{\mathbf{q}}^{0}$. It is immediate to see that if these satisfy (41) then they provide a feasible complementary slack solution of M-CLP/M-CLP∗.

(iii) This follows from Corollary S3.3, where a proof was given for $\lambda=\mu=0$. The result for $\lambda<0,\mu>0$ is essentially same.     ∎

Corollary 4.4**.**

Consider a solution of Internal-SCLP/SCLP∗ and of Boundary-LP/LP∗ that satisfies (41). Then for the solution of Boundary-LP/LP∗ the following holds:

(i) The solutions of Boundary-LP$(\tilde{U},\tilde{x})$ and of Boundary-LP${}^{*}(\tilde{P},\tilde{q})$ are complementary slack.

(ii) For $k\in{\cal{K}}^{=}$ we have $x_{k}^{\bullet}=0$ , and for $j\in{\cal{J}}^{=}$ we have $q_{j}^{\bullet}=0$.

(iii) For $j\in{\cal{J}}^{=}$ we have $\mathbf{u}_{j}^{N}=0$, and for $k\in{\cal{K}}^{=}$ we have $\mathbf{p}_{k}^{0}=0$.

Proof.

One can see that from (41) follows:

[TABLE]

In particular this implies:

(i) Variables $\mathbf{u}_{j}^{0},{\mathbf{q}}_{j}^{0}$, variables $\mathbf{p}_{k}^{N},{\mathbf{x}}_{k}^{N}$, variables $x_{k}^{\bullet},\mathbf{p}_{k}^{0}$ and variables $q_{j}^{\bullet},\mathbf{u}_{j}^{N}$ are complementary slack.

(ii) In the solution of Internal-SCLP/SCLP${}^{*}(x^{0},q^{N})$, $x_{k}^{0}=\tilde{x}_{k},\,k\in{\cal{K}}^{=}$, and hence $x_{k}^{\bullet}=0$ for $k\in{\cal{K}}^{=}$. Similarly, $q_{j}^{\bullet}=0$ for $j\in{\cal{J}}^{=}$.

(iii) By complementary slackness of $x^{0},\mathbf{p}^{0}$ and $q^{N},\mathbf{u}^{N}$, we also have that $\mathbf{p}_{k}^{0}=0$ for $k\in{\cal{K}}^{=}$, $\mathbf{u}_{j}^{N}=0$ for $j\in{\cal{J}}^{=}$.   ∎

The problems Boundary-LP(${\tilde{U}},{\tilde{x}}$), Boundary-LP${}^{*}({\tilde{P}},{\tilde{q}})$ may have non-unique optimal solutions for a given set of parameters ${\tilde{U}},{\tilde{x}},{\tilde{P}},{\tilde{q}}$. If there are multiple solutions then it is not obvious how to check if (41) holds for some solutions. The following Theorem answers these questions.

Theorem 4.5**.**

Under the Non-Degeneracy Assumption 1 there is a one to one correspondence between optimal solutions of M-CLP and Boundary-LP(${\tilde{U}},{\tilde{x}}$). The same holds for optimal solutions of M-CLP∗ and Boundary-LP${}^{*}({\tilde{P}},{\tilde{q}})$.

Proof.

From Proposition 4.2 it follows that under Non-Degeneracy Assumption 1 the formulation of Boundary-LP(${\tilde{U}},{\tilde{x}}$) problem is unique, whether M-CLP has a unique solution or not.

Furthermore, each optimal solution of M-CLP determines the boundary values $\mathbf{u}^{0},\mathbf{u}^{N},{\mathbf{x}}^{N}$ and the values of $x^{\bullet}=x^{0}-{\tilde{x}}$ that by Theorem 4.3(iii) is an optimal solution of Boundary-LP(${\tilde{U}},{\tilde{x}}$). Hence, each optimal solution of M-CLP uniquely defines an optimal solution of the Boundary-LP(${\tilde{U}},{\tilde{x}}$).

Conversely, each optimal solution of Boundary-LP(${\tilde{U}},{\tilde{x}}$) together with $\dot{x}(t),u(t)$ that are unique by Theorem D5.5(iii) produces a feasible solution of M-CLP. Moreover, objective values of these solutions satisfy:

[TABLE]

The second term on the right hand side is uniquely determined, and therefore it follows that all the optimal solutions of Boundary-LP(${\tilde{U}},{\tilde{x}}$) produces a feasible solution of M-CLP with same objective value. However, by Theorem 4.3 one of the optimal solutions of Boundary-LP(${\tilde{U}},{\tilde{x}}$), is obtained by the decomposition of an optimal solution of M-CLP and hence each optimal solution of Boundary-LP(${\tilde{U}},{\tilde{x}}$) uniquely defines an optimal solution of M-CLP.   ∎

Note: Theorem 4.5 explains why we choose the formulations (39), (40) of the Boundary-LP/LP∗ problems. This one to one correspondence of optimal solutions that is shown in the Theorem does not exist for the more natural pair of boundary problems (S13), (S14) in [31].

4.1 Boundary Simplex Dictionary

By the complementary slackness condition in equations (28), there are exactly $2K+2J$ boundary values which are 0, and up to $2K+2J$ positive ones. However, as it turns out, the number of positive primal boundary values may be different from $2K$. In fact the number of positive primal boundary values is $|{\cal{K}}_{0}|+|\overline{{\cal{J}}_{0}}|+|{\cal{K}}_{N+1}|+|\overline{{\cal{J}}_{N+1}}|=2K+L$, and correspondingly, the number of positive dual boundary values $|{\cal{J}}_{N+1}|+|\overline{{\cal{K}}_{N+1}}|+|{\cal{J}}_{0}|+|\overline{{\cal{K}}_{0}}|=2J-L$ where $L$ may be either 0 or positive or negative. In other words, the set of columns of \left[\begin{array}[]{cccc}A&0&I&0\\ A&A&0&I\end{array}\right] that correspond to positive primal boundary values may not be a basis, and similarly for the dual.

Recall also that the problems Boundary-LP$({\tilde{U}},{\tilde{x}})$ and Boundary-LP${}^{*}({\tilde{P}},{\tilde{q}})$ are not dual to each other. However, optimal solutions of Boundary-LP$({\tilde{U}},{\tilde{x}})$ and Boundary-LP${}^{*}({\tilde{P}},{\tilde{q}})$ that correspond to optimal solutions of M-CLP/M-CLP∗ are unique and are complementary slack. If possible we wish to find basic solutions that are summarized by a single boundary simplex dictionary. This dictionary will enables us to define a boundary pivot operation in Section 7. The following theorem states when this is possible.

Theorem 4.6**.**

If M-CLP/M-CLP∗ has a unique solution with $H_{\mathbb{P}}>0$, then its optimal base sequence satisfies $|{\cal{K}}^{\uparrow\downarrow}|+|\overline{{\cal{J}}_{0}}|+|{\cal{K}}_{N+1}|+|\overline{{\cal{J}}_{N+1}}|\leq 2K$ and $|{\cal{J}}_{0}|+|\overline{{\cal{K}}_{0}}|+|{\cal{J}}^{\uparrow\downarrow}|+|\overline{{\cal{K}}_{N+1}}|\leq 2J$.

Note: while $|{\cal{K}}_{0}|+|\overline{{\cal{J}}_{0}}|+|{\cal{K}}_{N+1}|+|\overline{{\cal{J}}_{N+1}}|$ may be larger than $2K$, the theorem states that under the given conditions, if $|{\cal{K}}_{0}|$ is replaced by $|{\cal{K}}^{\uparrow\downarrow}|$, the sum cannot exceed $2K$.

Proof.

If M-CLP/M-CLP∗ has a unique solution with $H_{\mathbb{P}}>0$, then the optimal solution of the corresponding Boundary-LP (39) has $|{\cal{K}}^{\uparrow\downarrow}|+|\overline{{\cal{J}}_{0}}|+|{\cal{K}}_{N+1}|+|\overline{{\cal{J}}_{N+1}}|$ positive variables. But if $|{\cal{K}}^{\uparrow\downarrow}|+|\overline{{\cal{J}}_{0}}|+|{\cal{K}}_{N+1}|+|\overline{{\cal{J}}_{N+1}}|>2K$ this implies that it has a non-basic optimal solution, and therefore its solution is non-unique. But then by Theorem 4.5 M-CLP has non-unique optimal solutions, which is a contradiction.   ∎

We now decompose the matrix of coefficients of (39) \left[\begin{array}[]{cccc}A&0&I&0\\ A&A&0&I\end{array}\right] into a basic and non-basic part, with similar decomposition for the matrix of coefficients of (40) \left[\begin{array}[]{cccc}A^{{\mbox{\tiny\bf\sf T}}}&0&-I&0\\ A^{{\mbox{\tiny\bf\sf T}}}&A^{{\mbox{\tiny\bf\sf T}}}&0&-I\end{array}\right]. Let ${\mathbf{B}}$ be a basic submatrix of the former, ${\mathbf{B}}^{*}$ a basic matrix of the latter, such that they are complementary slack, and let ${\mathbf{N}}$, ${\mathbf{N}}^{*}$ be the matrices composed from the corresponding non-basic columns, so that (as seen for example in [34]) ${\mathbf{B}}^{-1}{\mathbf{N}}=({\mathbf{B}}^{*-1}{\mathbf{N}}^{*})^{{\mbox{\tiny\bf\sf T}}}$. We then say that ${\mathbf{B}}$ and ${\mathbf{B}}^{*}$ are compatible.

We now consider a base sequence that satisfies $|{\cal{K}}^{\uparrow\downarrow}|+|\overline{{\cal{J}}_{0}}|+|{\cal{K}}_{N+1}|+|\overline{{\cal{J}}_{N+1}}|\leq 2K$ and $|{\cal{J}}^{\uparrow\downarrow}|+|\overline{{\cal{K}}_{N+1}}|+|{\cal{J}}_{0}|+|\overline{{\cal{K}}_{0}}|\leq 2J$. For such a base sequence we now specify rules to construct a pair of compatible bases of boundary variables:

define primal variables with indexes that are members of ${\cal{K}}^{\uparrow\downarrow},\overline{{\cal{J}}_{0}},{\cal{K}}_{N+1},$ $\overline{{\cal{J}}_{N+1}}$ as basic variables for ${\mathbf{B}}$ and dual variables with indexes that are members of ${\cal{J}}^{\uparrow\downarrow},\overline{{\cal{K}}_{N+1}},{\cal{J}}_{0},\overline{{\cal{K}}_{0}}$ as basic variables for ${\mathbf{B}}^{*}$,

• -

if the solution of M-CLP possesses $\geq 2K$ strictly positive primal boundary values then choose an arbitrary subset of size $|{\cal{K}}_{0}|+|\overline{{\cal{J}}_{0}}|+|{\cal{K}}_{N+1}|+|\overline{{\cal{J}}_{N+1}}|-2K$ of indexes out of ${\cal{K}}^{=}$ and define those $\mathbf{p}_{k}^{0}$ as dual basic variable in ${\mathbf{B}}^{*}$, and for the remaining indexes of ${\cal{K}}^{=}$ define $x_{k}^{\bullet}$ as primal basic variable in ${\mathbf{B}}$. Recall that for all the variables in ${\cal{K}}^{=}$, $x_{k}^{\bullet}=0$, and at the same time also $\mathbf{p}_{k}^{0}=0$ by complementary slackness. So there is no problem in naming either $x_{k}^{\bullet}$ as primal basic, or $\mathbf{p}_{k}^{0}$ as dual basic, for an arbitrary $k\in{\cal{K}}^{=}$.

Define all the variables $q_{j}^{\bullet},\,j\in{\cal{J}}^{=}$ as dual basic variables in ${\mathbf{B}}^{*}$.

• -

similarly, if the solution of M-CLP∗ possesses $\geq 2J$ strictly positive boundary values, choose an arbitrary subset of size $|{\cal{J}}_{N+1}|+|\overline{{\cal{K}}_{N+1}}|+|{\cal{J}}_{0}|+|\overline{{\cal{K}}_{0}}|-2J$ of indexes out of ${\cal{J}}^{=}$ and define those $\mathbf{u}_{j}^{N}$ as basic variables in ${\mathbf{B}}$, and for the remaining indexes of ${\cal{J}}^{=}$ define $q_{j}^{\bullet}$ as dual basic variable in ${\mathbf{B}}^{*}$. Define all the variables $x_{k}^{\bullet},\,j\in{\cal{K}}^{=}$ as primal basic variables in ${\mathbf{B}}$.

Using compatible boundary bases we define a boundary simplex dictionary ${\cal{D}}$ by:

[TABLE]

where ${\mathbf{v}}=\{x^{\bullet}=x^{0}-{\tilde{x}},\mathbf{u}^{0},{\mathbf{x}}^{N},\mathbf{u}^{N}\}$, ${\mathbf{v}}^{*}=\{q^{\bullet}=q^{N}-{\tilde{q}},\mathbf{p}^{N},{\mathbf{q}}^{0},\mathbf{p}^{0}\}$ are values of variables of Boundary-LP/LP∗ and ${\cal{B}},{\cal{B}}^{*}$ are sets of indexes of the corresponding basic variables.

In our algorithm, at iterations in which a boundary pivot will be required, changes in $({\cal{K}}_{n},{\cal{J}}_{n})_{n\in\{0,N+1\}}$ will be achieved by pivots of this specially constructed dictionary. In the proof of Theorem 7.1, in Appendix A.4 we show that the arbitrary choice involved in constructing the compatible bases for the dictionary does not lead to any ambiguity.

5 Validity Regions

Throughout this section we assume that Non-Degeneracy Assumption 1 holds. Similar to the definition of validity region S3.6 we define:

Definition 5.1**.**

Let $({\cal{K}}_{n},{\cal{J}}_{n})_{n=0}^{N+1}$ be a proper base sequence. Let ${\cal{V}}$ be the set of all $\rho=(\beta,\gamma,\lambda,\mu,T)$ for which this base sequence is optimal. Then ${\cal{V}}$ is called the validity region of $({\cal{K}}_{n},{\cal{J}}_{n})_{n=0}^{N+1}$. Let ${\cal{F}}$ be the union of all validity regions, then ${\cal{F}}$ is called the parametric feasible region.

As in Theorem S3.7, each validity region is a closed convex polyhedral cone (it may consist only of the origin). The parametric feasibility region is clearly convex, and it is a union of validity regions, hence ${\cal{F}}$ is also a closed convex polyhedral cone.

Similar to Corollary S3.8 we can state:

Corollary 5.2**.**

For a given optimal base sequence, within its validity region, the elements of $H_{\mathbb{P}}$ are affine functions of the non-zero elements of ${\mathbf{R}}$. Conversely, consider all optimal solutions $H$ that belong to a given optimal base sequence. Then the non-zero elements of ${\mathbf{R}}$ change linearly as a function of non-zero elements of $H_{\mathbb{P}}$.

Proposition 5.3**.**

The number of validity regions is bounded by $\binom{4(K+J)}{2(K+J)}\;2^{\binom{K+J}{K}}$.

Proof.

The total number of bases of the Rates-LP/LP∗ is $\binom{K+J}{K}$. Moreover, for a proper base sequence, under the non-degeneracy assumption, the objective values of the Rates-LP/LP∗ are strictly decreasing in $n$ (see Corollary 4 in [36]) and hence the number of sequences $({\cal{K}}_{n},{\cal{J}}_{n})_{n=1}^{N}$ is bounded by $2^{\binom{K+J}{K}}$. The number of $({\cal{K}}_{n},{\cal{J}}_{n})_{n\in\{0,N+1\}}$ is bounded by $\binom{4(K+J)}{2(K+J)}$ and hence the total number of proper base sequences is bounded by $\binom{4(K+J)}{2(K+J)}\;2^{\binom{K+J}{K}}$.   ∎

Note, a simpler somewhat larger bound is $2^{2^{K+J+1}}$.

Theorem 5.4**.**

For a proper base sequence $({\cal{K}}_{n},{\cal{J}}_{n})_{n=0}^{N+1}$ with validity region ${\cal{V}}$, and for boundary parameters $\rho$ in the interior of ${\cal{F}}$, the following statements are equivalent:

(i) $\rho$ is an interior point of ${\cal{V}}$,

(ii) $\rho$ does not belong to the validity region of any other proper base sequence,

(iii) the solution of (32) for the point $\rho$ is unique and satisfies $H_{\mathbb{P}}>0$.

Proof.

See Appendix A.3   ∎

Note that Theorem 5.4 implies that interiors of validity regions are disjoint.

Corollary 5.5**.**

If a proper base sequence has a validity region with non-empty interior, then the corresponding matrix $M$ is non-singular.

Proof.

By Theorem 5.4, at an interior point of the validity region the solution of (32) is unique, and hence $M$ is non-singular.   ∎

Corollary 5.6**.**

If a point $\rho$ is an interior point of some validity region then the solutions of M-CLP and M-CLP∗ for this point are unique.

Proof.

By Theorem 5.4 the sequence of bases is unique, $u(t),p(t),\dot{x}(t),\dot{q}(t)$ are unique by Theorem D5.5(iii), and so $M$ is unique, and the solution is unique.   ∎

6 Collisions

In this section we describe how the solution changes as the boundary parameters move from an internal point of a validity region to a boundary point of this validity region. We consider a validity region ${\cal{V}}$ with non-empty interior and base sequence $({\cal{K}}_{n},{\cal{J}}_{n})_{n=0}^{N+1}$. Let ${\cal{L}}(\theta)$ be a parametric line of boundary parameters, with ${\cal{L}}(\theta)$ in the interior of ${\cal{V}}$ for $\theta<\overline{\theta}$, ${\cal{L}}(\overline{\theta})$ on the boundary, with $H_{\mathbb{P}}(\theta)$ the solution at ${\cal{L}}(\theta)$. Then $H_{\mathbb{P}}(\theta)>0$ for $\theta<\overline{\theta}$, but some elements of $H_{\mathbb{P}}$ shrink to 0 as $\theta\nearrow\overline{\theta}$. We call this a collision.

We now list all types of collisions. For convenience we number them (a), (b), $\ldots$, (f) .

Internal collisions

(a)

State collision. A single value $x^{n}_{k},n=1,\ldots,N$, or $q^{n}_{j},n=0,\ldots,N-1$, shrinks to [math] (see Figure 9). This is classified in [36] as SCLP-III type collision.

(b)

Interval collision between non-adjacent bases. A sequence of one or more internal interval lengths, $\tau_{n_{*}},\ldots,\tau_{n_{**}}$ shrink to 0, where $1<n_{*}\leq n_{**}<N$ and $|B_{n_{*}-1}\setminus B_{n_{**}+1}|=2$, i.e. the basis preceding $n_{*}$ and the basis succeeding $n_{**}$ are not adjacent. This means that at the collision point on the boundary, at time $t=t_{n_{*}}=\cdots=t_{n_{**}}$ two different state variables hit the value 0 ($x_{k}(t)$ and $x_{k^{\prime}}(t)$, or $q_{j}(T-t)$ and $q_{j^{\prime}}(T-t)$ or $x_{k}(t)$ and $q_{j}(T-t)$ (see Figure 10)). This is classified in [36] as SCLP-II type collision.

(c)

Interval collision between adjacent bases or at 0 or $T$, variable does not become free. Includes the following three cases: (i) A sequence of one or more interval lengths, $\tau_{n_{*}},\ldots,\tau_{n_{**}}$ shrink to 0, where $1<n_{*}\leq n_{**}<N$ and $|B_{n_{*}-1}\setminus B_{n_{**}+1}|=1$, (ii) intervals $\tau_{1},\ldots,\tau_{n_{**}},$ $1\leq n_{**}<N$ shrink to zero, (iii) intervals $\tau_{n_{*}},\ldots,\tau_{N},\,1<n_{*}\leq N$ shrink to zero. In case (i), a state variable $x_{k}$ (or a dual state variable $q_{j}$) is positive in the intervals $t_{n_{*}-1}<t<t_{n_{*}}$ and $t_{n_{**}}<t<t_{n_{**}+1}$, but it is 0 at $t_{n_{*}}$ and at $t_{n_{**}}$. In the case (ii) the collision is at $t=0$ and a dual state variable $q_{j}$ in positive in $t_{n_{**}}<t<t_{n_{**}+1}$ but it is 0 at $t_{n_{**}}$. In case (iii) the collision is at $t=T$, and a state variable $x_{k}$ in positive in $t_{n_{*}-1}<t<t_{n_{*}}$ but it is 0 at $t_{n_{*}}$. These collisions are classified as SCLP-I type collisions in [36].

In addition, we also assume that in validity region ${\cal{V}}$, for the state variable $x_{k}(t)$ (or $q_{j}(t)$) there is at least one time interval in $[0,t_{n_{*}}]$ or in $[t_{n_{**}},T]$ such that $x_{k}(t)=0$ (or $q_{j}(T-t)=0$) in that interval. This implies that the variable $x_{k}(t)$ (or $q_{j}(t)$) would be tied even if $x_{k}(t)$ (or $q_{j}(T-t)$) would be positive in the whole interval $[t_{n_{*}},t_{n_{**}}]$. This additional condition is what is meant by saying the variable involved in the collision does not become free. In Figure 11 a single interval $\tau_{n^{*}}=\tau_{n^{**}}=\tau_{4}$ shrinks to zero, but the variable $x_{k}(t)$ remains tied, since it satisfies $x_{k}(t)=0,\,0<t<t_{1}$.

Collisions that involve boundary values

(d)

Interval collision between adjacent bases or at 0 or $T$, variable become free. This is the same as (c) except that the variable $x_{k}$ (or $q_{j}$) become free in the sense that in the interior of the validity region ${\cal{V}}$, in the time intervals $(0,t_{n_{*}})$ and $(t_{n_{**}},T)$, $x_{k}(t)>0$ (or $q_{j}(T-t)>0$) (see Figure 12).

(e)

Boundary collision. One of the $>0$ boundary values including any of $\mathbf{u}^{0}_{j},\mathbf{u}^{N}_{j},\mathbf{p}^{0}_{k},\mathbf{p}^{N}_{k},{\mathbf{x}}_{k}^{N},$ ${\mathbf{q}}_{j}^{0}$ shrinks to zero, or one of the values $x_{k}^{0},\,k\in{\cal{K}}^{\uparrow\downarrow}$ or $q_{j}^{N},\,j\in{\cal{J}}^{\uparrow\downarrow}$, shrinks to 0, or one of the values of $x_{k}^{0},\,k\in{\cal{K}}^{=}$ for which $\dot{x}_{k}^{1}>0$ (see Figure 13) or one of the values of $q_{j}^{N},\,j\in{\cal{J}}^{=}$ for which $\dot{q}_{j}^{N}>0$ shrinks to 0.

(f)

Joint collision. One of the values of $x_{k}^{0},\,k\in{\cal{K}}^{=}$ for which $\dot{x}_{k}^{1}<0$ shrinks to 0, and as a result the intervals $\tau_{1},\ldots,\tau_{n_{*}}$ shrink to zero, where in the interior of the validity region $x_{k}(t)>0$ for $0<t<t_{n_{*}}$, and $x_{k}(t_{n_{*}})=0$ (see Figure 14).

Alternatively, one of the values of $q_{j}^{N},\,j\in{\cal{J}}^{=}$ for which $\dot{q}_{j}^{N}<0$ shrinks to 0, and as a result the intervals $\tau_{n_{**}},\ldots,\tau_{N}$ shrink to zero, where in the interior of the validity region $q_{j}(T-t)>0$ for $t_{n_{**}-1}<t<T$, and $q_{j}(T-t_{n_{**}-1})=0$.

Multiple collisions

Any collision that is not one of the types (a)–(f) is a multiple collision. Such collisions consists of a combination of several collisions of types (a)–(f). We will exclude the possibility of multiple collisions in Sections 7–9, but will return to discuss them in Appendix A.1.

7 M-CLP pivots

An M-CLP pivot uses the solution in validity region ${\cal{V}}$ with non-empty interior to construct the solution in a neighboring validity region with non-empty interior ${\cal{W}}$. Let ${\cal{L}}(\theta),\,0<\theta<1$ be a parametric line of boundary parameters, such that ${\cal{L}}(\theta)\in\mbox{interior of }{\cal{V}},\,\theta<\overline{\theta}$, ${\cal{L}}(\theta)\in\mbox{interior of }{\cal{W}},\,\theta>\overline{\theta}$, and ${\cal{L}}(\overline{\theta})$ is on the boundary of ${\cal{V}}$ and of ${\cal{W}}$. We assume that at the collision point ${\cal{L}}(\overline{\theta})$, the collision from ${\cal{V}}$ as well as the collision from ${\cal{W}}$ are both of types (a)–(f). In an M-CLP pivot we will start from the region ${\cal{V}}$ with optimal base sequence $({\cal{K}}_{n},{\cal{J}}_{n})_{n=0}^{N+1}$, and with a specified type of collision, and we will show how to retrieve the collision from the other side, and from it construct the base sequence $({\cal{K}}^{\prime}_{n},{\cal{J}}^{\prime}_{n})_{n=0}^{N^{\prime}+1}$ optimal in ${\cal{W}}$.

There are three types of pivots:

Pivots that involve only changes from $({\cal{K}}_{n},{\cal{J}}_{n})_{n=1}^{N}$ to $({\cal{K}}^{\prime}_{n},{\cal{J}}^{\prime}_{n})_{n=1}^{N^{\prime}}$, which are pivots on the Internal-SCLP, similar to SCLP pivots described in [36]. We refer to these as internal pivots.

• -

Pivots that involve changes from $({\cal{K}}_{n},{\cal{J}}_{n})_{n\in\{0,N+1\}}$ to $({\cal{K}}^{\prime}_{n},{\cal{J}}^{\prime}_{n})_{n\in\{0,N^{\prime}+1\}}$ in which the solution at the collision point is unique. We refer to these as boundary pivots of type I.

• -

Pivots that involve changes from $({\cal{K}}_{n},{\cal{J}}_{n})_{n\in\{0,N+1\}}$ to $({\cal{K}}^{\prime}_{n},{\cal{J}}^{\prime}_{n})_{n\in\{0,N^{\prime}+1\}}$ in which the solution at the collision point is not unique. We refer to these as boundary pivots of type II.

We now describe all the possible types of pivots:

Pure Internal-SCLP pivots

When the collision from validity region ${\cal{V}}$ is of types (a), (b), or (c) the pivot is an internal pivot.

Collision of type (a): A single state variable $x^{n}_{k}$, or $q^{n}_{j}$ shrinks to 0. This is classified as SCLP-III type collision in [36]. In this pivot new internal bases $D_{1},\ldots,D_{M}$ are inserted before $B_{1}$ (if $n=0$), or after $B_{N}$ (if $n=N$), or between $B_{n}$ and $B_{n+1}$ (otherwise). At the collision point these new bases have interval length 0, and they become $>0$ in the interior of ${\cal{W}}$. The collision from ${\cal{W}}$ is of type (c) if the shrinking state variable is tied, or of type (d) if it is free.

• -

Collision of type (b): Intervals $\tau_{n^{*}},\ldots,\tau_{n^{**}},\,1<n^{*}\leq n^{**}<N$ shrink to zero, with $|B_{n_{*}-1}\setminus B_{n_{**}+1}|=2$. This is classified as SCLP-II type collision in [36]. The internal pivot involves removing the bases $B_{n_{*}},\ldots,B_{n_{**}}$, and inserting new bases $D_{1},\ldots,D_{M},\,M\geq 1$ in their place. At the collision point these new bases have interval length 0, and they become $>0$ in the interior of ${\cal{W}}$. The collision from ${\cal{W}}$ is also of type (b).

• -

Collision of type (c): Intervals $\tau_{n^{*}},\ldots,\tau_{n^{**}}$ shrink to zero, where either $1<n^{*}\leq n^{**}<N$ with $|B_{n_{*}-1}\setminus B_{n_{**}+1}|=1$, or else $n^{*}=1$, or else $n^{**}=N$, but no tied state variable becomes free. This is classified as SCLP-I type collision in [36]. The internal pivot involves removing the bases $B_{n_{*}},\ldots,B_{n_{**}}$. The collision from ${\cal{W}}$ is of type (a).

M-CLP pivots that involve boundary pivots

This includes collisions of type (d), (e), (f) from the ${\cal{V}}$ side, and the pivot consists of up to three steps: it may require an initial Internal-SCLP pivot (pre-boundary step), it may then require a pivot on the boundary dictionary (boundary step), and it may finally require another Internal-SCLP pivot (post-boundary step).

The pre-boundary step is determined by the type of collision on the ${\cal{V}}$ side. Once pre-boundary SCLP pivot is complete the boundary pivot is carried out as follows: First the Boundary-LP simplex dictionary for ${\cal{V}}$ is constructed. Then the leaving variable is determined according to the type (d), (e) or (f). It is then determined if the pivot is type I or type II. If it is type I no pivot is performed, but the bases $({\cal{K}}_{n},{\cal{J}}_{n})_{n\in\{0,N+1\}}$ may change. If it is type II, a pivot of the boundary dictionary is performed, and the new bases $({\cal{K}}^{\prime}_{n},{\cal{J}}^{\prime}_{n})_{n\in\{0,N^{\prime}+1\}}$ are obtained. This also determines the type of collision from the ${\cal{W}}$ side. Finally, according to the type of collision from the ${\cal{W}}$ side a post-boundary SCLP pivot may be performed. We now describe these steps.

The pre-boundary step

Collision of type (d) or (f): The bases of intervals that shrink to 0 are removed.

• -

Collision of type (e): No pre-boundary pivot is required.

We denote by $({\tilde{{\cal{K}}}}_{n},{\tilde{{\cal{J}}}}_{n})_{n=1}^{N^{\prime\prime}}$ the base sequence obtained from $({\cal{K}}_{n},{\cal{J}}_{n})_{n=1}^{N}$ after the pre-boundary step, and refer to it as the intermediate base sequence.

Constructing the Boundary-LP simplex dictionary for ${\cal{V}}$:

In all three cases we follow the rules detailed in Section 4.1. In case (d) we need some additional rules: If $x_{k}$ became free then if $|{\cal{J}}^{\uparrow\downarrow}|+|\overline{{\cal{K}}_{N+1}}|+|{\cal{J}}_{0}|+|\overline{{\cal{K}}_{0}}|<2J$ we choose $\mathbf{p}_{k}^{0}$ as a dual basic in ${\mathbf{B}}^{*}$, and if $|{\cal{J}}^{\uparrow\downarrow}|+|\overline{{\cal{K}}_{N+1}}|+|{\cal{J}}_{0}|+|\overline{{\cal{K}}_{0}}|=2J$, we choose $x_{k}^{\bullet}$ as basic in ${\mathbf{B}}$. Similarly, if $q_{j}$ becomes free then if $|{\cal{K}}^{\uparrow\downarrow}|+|\overline{{\cal{J}}_{0}}|+|{\cal{K}}_{N+1}|+|\overline{{\cal{J}}_{N+1}}|<2K$ we choose $\mathbf{u}_{j}^{N}$ as basic in ${\mathbf{B}}$, and if $|{\cal{K}}^{\uparrow\downarrow}|+|\overline{{\cal{J}}_{0}}|+|{\cal{K}}_{N+1}|+|\overline{{\cal{J}}_{N+1}}|=2K$ we choose $q_{j}^{\bullet}$ as dual basic in ${\mathbf{B}}^{*}$. We denote the Boundary-LP simplex dictionary for ${\cal{V}}$ by ${\cal{D}}$.

Determining the leaving variable

In cases (e), (f) the variable leaving the basis is the boundary variable that has shrunk to 0 in the collision from ${\cal{V}}$. We denote it by ${\mathbf{v}}$, and it can be either a primal variable or a dual variable. We denote its dual variable by ${\mathbf{v}}^{*}$.

In case (d), if $x_{k}$ has become free and $x_{k}^{\bullet}$ is basic, or if $q_{j}$ has become free and $q_{j}^{\bullet}$ is basic, then no variable leaves the basis and there is no pivot. We classify this case as a pivot of type I. Otherwise ${\mathbf{v}}=\mathbf{p}_{k}^{0}$ or ${\mathbf{v}}=\mathbf{u}_{j}^{N}$ respectively.

Determine the type of boundary pivot and the entering variable

We have the following possibilities:

If ${\mathbf{v}}$ is a primal (or a dual) variable, examine the dual (or the primal) variables which are ${\mathbf{w}}^{*}=0$. If any of them has a non-zero pivot element (i.e. if ${\mathbf{v}}$ is in row $i$ and ${\mathbf{w}}^{*}$ is in column $j$, the element ${\hat{{\mathbf{A}}}}_{ij}$) then this is a type I pivot.

• -

Otherwise the entering variable ${\mathbf{w}}$ is determined, as in standard LP, by the ratio test (this is the ratio test of the simplex method of linear programming, to determine which variable becomes tight and leaves the primal or the dual basis).

• –

If there is a single candidate variable ${\mathbf{w}}$ to enter with ratio $>0$ then we have a pivot of type II and ${\mathbf{w}}$ is the entering variable.

• –

If there are multiple candidates to enter with ratio $>0$, then there are several variables shrinking to [math] on the ${\cal{W}}$ side. This indicates that there is a multiple collision from the ${\cal{W}}$ side. We discuss multiple collisions in Appendix A.1.

Boundary-pivot of type I:

The solution at the collision point is unique, and there is no need for a pivot on the dictionary. In the case of collision type (d), when ${\mathbf{v}}=\mathbf{p}_{k}^{0}$ (or ${\mathbf{v}}=\mathbf{u}_{j}^{N}$) we have $({\cal{K}}_{n},{\cal{J}}_{n})_{n\in\{0,N+1\}}=({\cal{K}}^{\prime}_{n},{\cal{J}}^{\prime}_{n})_{n\in\{0,N^{\prime}+1\}}$, and the only change is that $x_{k}$ (or $q_{j}$) is no longer tied but becomes free in ${\cal{W}}$. Otherwise (in case of collision type (e) or (f)), the solution at the collision point is unique, but $({\cal{K}}_{n},{\cal{J}}_{n})_{n\in\{0,N+1\}}$ $\neq({\cal{K}}^{\prime}_{n},{\cal{J}}^{\prime}_{n})_{n\in\{0,N^{\prime}+1\}}$.

Boundary-pivot of type II:

Assume ${\mathbf{v}}$ leaving and ${\mathbf{w}}$ entering. Let ${\mathbf{v}}^{*}$, ${\mathbf{w}}^{*}$ be their dual variables. Then the effect of the pivot is as follows: Because the leaving variable ${\mathbf{v}}=0$, after the pivot the new basic variable ${\mathbf{w}}=0$. However, since the test ratio is $>0$, the variable ${\mathbf{w}}^{*}$ that leaves the dual basis is $>0$, and when it leaves the basis its value will jump down to 0. Also, the variable entering the dual basis ${\mathbf{v}}^{*}$ will jump from the non-basic value of 0 to a value $>0$. In addition, all the basic boundary variables in the same category as ${\mathbf{w}}^{*},{\mathbf{v}}^{*}$ (that is in ${\mathbf{B}}^{*}$ or in ${\mathbf{B}}$ according to whether ${\mathbf{v}}$ is primal or dual) will change their values.

Determine the variable shrinking to [math] from the ${\cal{W}}$ side:

From the ${\cal{W}}$ side we may have a single boundary variable shrinking to [math], which we denote by ${\mathbf{v}}^{\prime}$, or there may be no boundary value that shrinks to 0 and we then define ${\mathbf{v}}^{\prime}=\emptyset$.

If the pivot is type I and the collision from ${\cal{V}}$ is (d) then ${\mathbf{v}}^{\prime}=\emptyset$.

• -

If the pivot is type I and the collision from ${\cal{V}}$ is (e) or (f) then ${\mathbf{v}}^{\prime}={\mathbf{v}}^{*}$.

• -

In case of boundary pivot of type II, if ${\mathbf{w}}^{*}=x^{\bullet}_{k^{\prime}},{\tilde{x}}_{k^{\prime}}>0$ (or if ${\mathbf{w}}^{*}=q^{\bullet}_{j^{\prime}},{\tilde{q}}_{j^{\prime}}>0$) then ${\mathbf{v}}^{\prime}=\emptyset$.

• -

Otherwise, for pivot of type II, ${\mathbf{v}}^{\prime}={\mathbf{w}}$.

Determine the new boundary basis:

To determine the new boundary bases $({\cal{K}}^{\prime}_{n},$ ${\cal{J}}^{\prime}_{n})_{n\in\{0,N^{\prime}+1\}}$ we first find the new values of $x^{0},q^{N}$, from the new values of $x^{\bullet}$, $q^{\bullet}$ after the pivot, using $x^{0}=x^{\bullet}+{\tilde{x}},q^{N}=q^{\bullet}+{\tilde{q}}$. Recall that ${\tilde{x}},\,{\tilde{q}}$ are uniquely determined at ${\cal{L}}(\overline{\theta})$ by Proposition 4.2.

Following this step we define $({\cal{K}}^{\prime}_{n},{\cal{J}}^{\prime}_{n})_{n\in\{0,N^{\prime}+1\}}$ as follows:

[TABLE]

This simply includes all $>0$ boundary values, and the variable ${\mathbf{v}}^{\prime}$ (if ${\mathbf{v}}^{\prime}\neq\emptyset$), and includes exactly $2(K+J)$ elements.

Test for a multiple collision from the ${\cal{W}}$ side:

If ${\mathbf{w}}^{*}=x_{k}^{\bullet}$ or ${\mathbf{w}}^{*}=q_{j}^{\bullet}$ and at ${\cal{L}}(\overline{\theta})$ the minimum of $x_{k}(t)$ or $q_{j}(t)$ is obtained in more than one point, as seen from $|\arg\min_{0\leq n\leq N}x^{n}_{k}|>1$ (or from $|\arg\min_{0\leq n\leq N}q^{n}_{j}|>1$), then there is a multiple collision from the ${\cal{W}}$ side. This is illustrated in Figure 15 where $x_{k}^{\bullet}(\theta)$ denotes the value of the optimal $x_{k}^{\bullet}$ for parameters ${\cal{L}}(\theta)$.

Determine the type of collision from the ${\cal{W}}$ side:

We assume that no multiple collision is discovered in the previous step.

If ${\mathbf{v}}^{\prime}=\emptyset$ and the boundary pivot is of type I, the collision from the ${\cal{W}}$ side is type (a).

• -

If ${\mathbf{v}}^{\prime}=\emptyset$ and the boundary pivot is of type II, the collision form the ${\cal{W}}$ side is type (d).

• -

If ${\mathbf{v}}^{\prime}=x^{\bullet}_{k^{\prime}}$ and ${\cal{K}}^{\prime}_{0}$ is not compatible with ${\tilde{{\cal{K}}}}_{1}$ (of the intermediate base sequence), or if ${\mathbf{v}}^{\prime}=q^{\bullet}_{j^{\prime}}$ and ${\cal{J}}^{\prime}_{N+1}$ is not compatible with ${\tilde{{\cal{J}}}}_{N}$, then the collision from the ${\cal{W}}$ side is of type (f).

• -

For any other ${\mathbf{v}}^{\prime}\neq\emptyset$, the collision from the ${\cal{W}}$ side is of type (e)

The post-boundary step

In the post-boundary step we perform an SCLP pivot and build a base sequence $({\cal{K}}^{\prime}_{n},{\cal{J}}^{\prime}_{n})_{n=1}^{N^{\prime}}$ that is optimal on the ${\cal{W}}$ side.

If the collision from the ${\cal{W}}$ side is of type (d), and ${\mathbf{w}}^{*}=x_{k}^{\bullet}$ and $t_{n}$ is the minimum point of $x_{k}(t)$ at ${\cal{L}}(\overline{\theta})$, or ${\mathbf{w}}^{*}=q_{j}^{\bullet}$ and $t_{n}$ is the minimum point of $q_{j}(T-t)$ at ${\cal{L}}(\overline{\theta})$, the bases $D_{1},\ldots,D_{M}$ need to be inserted between $\tilde{B}_{n},\tilde{B}_{n+1}$, in the sequence $({\tilde{{\cal{K}}}}_{n},{\tilde{{\cal{J}}}}_{n})_{n=1}^{N^{\prime\prime}}$ as in SCLP-III type collision in [36].

• -

If the collision from the ${\cal{W}}$ side is of type (f), then if ${\mathbf{v}}^{\prime}=x_{k}^{\bullet}$, then a number of bases $D_{1},\ldots,D_{M}$ are inserted before $\tilde{B}_{1}$, and if ${\mathbf{v}}^{\prime}=q_{j}^{\bullet}$, then a number of bases $D_{1},\ldots,D_{M}$ are inserted after $\tilde{B}_{N^{\prime\prime}}$, as in SCLP-IV type collision in [36].

We now formulate the following two theorems that state the conditions under which the pivot operation is well defined and constructs an optimal solution for the neighboring region ${\cal{W}}$.

Theorem 7.1**.**

If the collisions on both the ${\cal{V}}$ and the ${\cal{W}}$ side are of types (a)–(f) then the M-CLP pivot constructs a unique base sequence $({\cal{K}}^{\prime}_{n},{\cal{J}}^{\prime}_{n})_{n=0}^{N^{\prime}+1}$.

Proof.

See Appendix A.4   ∎

Theorem 7.2**.**

If the collisions on both the ${\cal{V}}$ side as $\theta\nearrow\overline{\theta}$ and on the ${\cal{W}}$ side as $\theta\searrow\overline{\theta}$ are both of types (a)–(f) then the base sequence $({\cal{K}}^{\prime}_{n},{\cal{J}}^{\prime}_{n})_{n=0}^{N^{\prime}+1}$ constructed by the M-CLP pivot is the optimal base sequence in ${\cal{W}}$.

Proof.

See Appendix A.4   ∎

8 Algorithm under simplifying assumption

Consider an M-CLP problem and let $\beta_{g},\gamma_{g},$ $\lambda_{g},\mu_{g},T_{g}$ be a goal set of boundary parameters in the feasible parametric region, i.e. the set of boundary parameters for which an optimal solution should be found. Let $\beta_{0},\gamma_{0},\lambda_{0},\mu_{0},T_{0}$ be another set of boundary parameters with a known optimal solution. Define the parametric line of boundary parameters as ${\cal{L}}(\theta)=(1-\theta)(\beta_{0},\gamma_{0},\lambda_{0},\mu_{0},T_{0})+\theta(\beta_{g},\gamma_{g},\lambda_{g},\mu_{g},T_{g})$, $0\leq\theta\leq 1$. The algorithm described in this section requires the following assumption:

Assumption 2** (Simplifying Assumption).**

Assume that all collisions along the parametric line ${\cal{L}}(\theta)$ are of type (a)-(f).

We discuss the algorithm without this assumption in Appendix A.1

8.1 Properties of the parametric line

Theorem 8.1**.**

Let ${\cal{L}}(\theta),0\leq\theta\leq 1$ be a parametric line of boundary parameters. If M-CLP/M-CLP∗ are feasible for both goal set of parameters ${\cal{L}}(1)$ and initial set of parameters ${\cal{L}}(0)$, then:

(i) all points of the parametric line belong to the parametric-feasible region,

(ii) if ${\cal{L}}(0)$ is an interior point of some validity region, then all points of the parametric line except possibly ${\cal{L}}(1)$ belong to the interior of the parametric-feasible region,

(iii) the parametric line can be partitioned to a minimal finite set of closed intervals of positive length $[\theta_{\ell-1},\theta_{\ell}],$ $\ell=1,\dots,L\leq\binom{4(K+J)}{2(K+J)}\;2^{\binom{K+J}{K}}$ such that for each interval there is a (possibly non-unique) proper base sequence which is optimal for every ${\cal{L}}(\theta),\,\theta\in[\theta_{\ell-1},\theta_{\ell}]$.

Proof.

See Appendix A.4   ∎

Proposition 8.2**.**

Consider a continuous parametric line of boundary parameters ${\cal{L}}(\theta)$. Then the decomposition parameters ${\tilde{U}}(\theta),{\tilde{x}}(\theta),{\tilde{P}}(\theta),{\tilde{q}}(\theta)$ are continuous piecewise affine functions of $\theta$.

Proof.

Within each closed validity region ${\tilde{U}}(\theta),{\tilde{x}}(\theta),{\tilde{P}}(\theta),{\tilde{q}}(\theta)$ are affine by Corollary 5.2. Moreover, by Proposition 4.2, for every set of boundary parameters ${\tilde{U}},{\tilde{x}},{\tilde{P}},{\tilde{q}}$ are unique, and therefore ${\tilde{U}}(\theta),{\tilde{x}}(\theta),{\tilde{P}}(\theta),{\tilde{q}}(\theta)$ are continuous at the boundary points of the validity regions ${\cal{L}}(\theta_{\ell})$.     ∎

8.1.1 Initial solution

The following theorem shows how to obtain a set of boundary parameters for a simple initial optimal solution.

Theorem 8.3**.**

Consider a set of boundary parameters $\beta_{0},\gamma_{0},\lambda_{0},\mu_{0},T_{0}$ that satisfy:

[TABLE]

then M-CLP/M-CLP∗ have a single-interval optimal solution, with $U(t)=P(t)=0$, and $x(t),-q(t)$ equal to the right hand sides of (3), (3), and this set of boundary values is an interior point of the validity region.

Proof.

It is easy to see that under (44) $U(t)=0,\,x(t)=\beta_{0}+bt,\,0\leq t<T_{0},\,x(T)=\beta_{0}+bT_{0}+\lambda_{0},\,P(t)=0,\,q(t)=-\gamma_{0}-ct,0\leq t<T_{0},$ $q(T)=-\gamma_{0}-cT_{0}-\mu_{0}$ are feasible complementary slack solutions of M-CLP/M-CLP∗, and hence are optimal. They define an optimal base sequence $({\cal{K}}_{n},{\cal{J}}_{n})_{n=0}^{2},\,{\cal{K}}_{0}={\cal{K}}_{1}={\cal{K}}_{2}=\{1,\dots,K\},\,{\cal{J}}_{0}={\cal{J}}_{1}={\cal{J}}_{2}=\{1,\dots,J\}$, with a single interval solution. Moreover, one can see that the solution to the corresponding Boundary-LP/LP∗ is unique, and hence by Theorem 4.5 this solution of M-CLP/M-CLP∗ is unique. It is easy to see that $H_{\mathbb{P}}>0$ and so by Theorem 5.4 the parameter values are an interior point of its validity region.     ∎

8.2 The algorithm

We now describe a simplex-type algorithm to solve M-CLP/M-CLP∗ in a finite number of steps. For the algorithm to work we need the Non-Degeneracy Assumption 1 and the Simplifying Assumption 2 to hold.

8.2.1 Input

Problem data consists of $A,b,c$, and goal parameters ${\cal{L}}(1)=(\beta_{g},\gamma_{g},\lambda_{g}\leq 0,\mu_{g}\geq 0,T_{g})$.

8.2.2 Output

Optimal solution of M-CLP/M-CLP∗ at ${\cal{L}}(1)$: optimal base sequence, breakpoint times, impulse controls at [math] and $T_{g}$, rates $u,\dot{x},p,\dot{q}$ in each interval, states $x(t),q(T_{g}-t)$ at each breakpoint.

8.2.3 Feasibility test

Solve Test-LP/LP∗ (23) for ${\cal{L}}(1)$ to determine whether both M-CLP/M-CLP∗ are feasible, or both are infeasible, or one is infeasible and the other unbounded. Proceed to solve the problems if both are feasible, otherwise stop.

8.2.4 Initialization

Choose ${\cal{L}}(0)=(\beta_{0},\gamma_{0},\lambda_{0},\mu_{0},T_{0})$ that satisfy:

[TABLE]

Define $\delta{\cal{L}}={\cal{L}}(1)-{\cal{L}}(0)$.

The initial solution has $U(t)=P(t)=0$, $x(t),-q(t)$ are r.h.s. of (3), (3) at ${\cal{L}}(0)$, and the initial optimal base sequence is $({\cal{K}}_{n},{\cal{J}}_{n})_{n=0}^{2}=\{({\cal{K}}_{n}=\{1,\dots,K\},{\cal{J}}_{n}=\{1,\dots,J\}),n=0,1,2\}$.

Set: $\ell:=1$,  $\theta_{\ell}:=0$,  $({\cal{K}},{\cal{J}})^{(\ell)}:=({\cal{K}}_{n},{\cal{J}}_{n})_{n=0}^{2}$.

8.2.5 Iteration

$\bullet$ Set up equations

Calculate rates for all bases in $({\cal{K}},{\cal{J}})^{(\ell)}$, construct the coefficients of the matrix $M_{\ell}$ for $({\cal{K}},{\cal{J}})^{(\ell)}$. Construct ${\mathbf{R}}_{\ell}$, $\delta{\mathbf{R}}_{\ell}$ from ${\cal{L}}(\theta_{\ell})$ and $\delta{\cal{L}}$.

$\bullet$ Calculate current solution and its $\theta$ gradient

Solve the equations

[TABLE]

$\bullet$ Find right endpoint of the validity region

Find

[TABLE]

The right endpoint of the validity region is $\overline{\theta}:=\theta_{\ell}+\overline{\Delta}$.

$\bullet$ Termination test and conclusion

If $\theta_{\ell}<1\leq\overline{\theta}$ terminate: $({\cal{K}},{\cal{J}})^{(\ell)}$ is the optimal base-sequence, current rates are optimal, remaining output is obtained from $H=H^{\ell}+(1-\theta_{\ell})\delta H^{\ell}$.

$\bullet$ Update solution

Calculate values of the solution at $\overline{\theta}$, as $\overline{H}=H^{\ell}+\overline{\Delta}\delta H^{\ell}$.

$\bullet$ Classify collision as $\theta\nearrow\overline{\theta}$

Find the set of elements of $H^{\ell}$ that are 0 at ${\cal{L}}(\overline{\theta})$. Split ${\cal{K}}_{0},{\cal{J}}_{N+1}$ to sets ${\cal{K}}^{=},{\cal{K}}^{\uparrow\downarrow},{\cal{J}}^{=},{\cal{J}}^{\uparrow\downarrow}$. Classify the collision to types (a)-(f) as discussed in Section 6.

$\bullet$ Pivot

Perform M-CLP pivot to determine a new base sequence, $({\cal{K}}^{\prime}_{n},{\cal{J}}^{\prime}_{n})_{n=0}^{N^{\prime}+1}$ as discussed in Section 7. Set $\ell:=\ell+1$, $\theta_{\ell}:=\overline{\theta}$, $({\cal{K}},{\cal{J}})^{(\ell)}:=({\cal{K}}^{\prime}_{n},{\cal{J}}^{\prime}_{n})_{n=0}^{N^{\prime}+1}$ and move to the next iteration.

8.3 Verification of the algorithm

In this section we verify that under Simplifying Assumption 2 and Non-Degeneracy Assumption 1 all calculations described in Section 8 are possible, and for any pair of feasible M-CLP/M-CLP∗ problems the algorithm produces an optimal solution in a finite bounded number of iterations. The main point here is to show that the new base sequence is optimal for $\theta>\overline{\theta}$ small enough, and we show that this property follows from Theorem 7.2.

Theorem 8.4**.**

Assume that M-CLP/M-CLP∗ are feasible for the goal set of parameters $\beta_{g},\gamma_{g},T_{g},$ $\lambda_{g},\mu_{g}$. Assume also Non-Degeneracy Assumption 1 and Simplifying Assumption 2. Then the algorithm will compute the solution of M-CLP/M-CLP${}^{*}(\beta_{g},\gamma_{g},T_{g},\lambda_{g},\mu_{g})$ in a finite number of iterations, bounded by $L\leq\binom{4(K+J)}{2(K+J)}\;2^{\binom{K+J}{K}}$.

Proof.

The proof follows from the following sequence of Propositions.   ∎

Proposition 8.5**.**

The feasibility test can be performed, to determine whether M-CLP is feasible and bounded, or is infeasible, or is unbounded.

Proof.

Clearly, one can determine if Test-LP (23) is feasible and bounded, or is infeasible or is unbounded, which determines the same for M-CLP.     ∎

Proposition 8.6**.**

Initialization phase of the algorithm yields a unique optimal solution with a unique optimal base sequence, and ${\cal{L}}(0)$ is an interior point of its validity region.

Proof.

Follows from Theorem 8.3.     ∎

Proposition 8.7**.**

Assume that ${\cal{L}}(\theta)$ for some $\theta>\theta_{\ell}$ is in the interior of the validity region of $({\cal{K}},{\cal{J}})^{(\ell)}$. Then setting up equations, calculation of the current solution and its gradient, and calculation of the right endpoint $\overline{\theta}$, as well as the termination test can be carried out.

Proof.

The base sequence $({\cal{K}},{\cal{J}})^{(\ell)}$ determines the structure of $M_{\ell}$ and of ${\mathbf{R}}_{\ell}$. Solution of the Rates-LP/LP∗ for $({\cal{K}}_{n},{\cal{J}}_{n})_{n=1}^{N}$ is unique by non-degeneracy assumption, and provides the non-zero coefficients of $M_{\ell}$. The non-zero coefficients of ${\mathbf{R}}_{\ell}$ are given by ${\cal{L}}(\theta_{\ell})$. $\delta{\mathbf{R}}_{\ell}$ has the same structure as ${\mathbf{R}}_{\ell}$, and its non-zero elements are given by $\delta{\cal{L}}$.

We assume in the statement of the proposition that the base sequence $({\cal{K}},{\cal{J}})^{(\ell)}$ has non-empty interior, and hence by Proposition 5.5 the matrix $M_{\ell}$ is non-singular. Hence $H^{\ell}$ is the unique solution of $M_{\ell}H^{\ell}={\mathbf{R}}_{\ell}$ and $\delta H^{\ell}$ is the unique solution of $M_{\ell}\delta H^{\ell}=\delta{\mathbf{R}}_{\ell}$. We denote by $H(\theta)$ the solution at point ${\cal{L}}(\theta)$. Then at the point $\theta>\theta_{\ell}$ in the interior of the validity region the solution is determined by

[TABLE]

That $\overline{\Delta}>0$ follows because the interior of the validity region of $({\cal{K}},{\cal{J}})^{(\ell)}$ is non-empty and contains an interval of the line ${\cal{L}}(\theta)$. It is possible to have $1/\overline{\Delta}=0$ and $\overline{\Delta}=\infty$.

Finally, if any of $\delta H_{r}<0$ then by $H_{\mathbb{P}}(\theta)\geq 0$ it follows that $H_{r}>0$, and the ratio $-\delta H_{r}/H_{r}>0$ and it is finite, and so $1/\overline{\Delta}>0$ and finite so that $0<\overline{\Delta}<\infty$. Otherwise $1/\overline{\Delta}=0$ and $\overline{\Delta}=\infty$.     ∎

Proposition 8.8**.**

M-CLP pivot can be performed and is unique. There exists $\theta>\overline{\theta}$ small enough, such that base sequence $({\cal{K}},{\cal{J}})^{(\ell+1)}$ is optimal for this $\theta$ and ${\cal{L}}(\theta)$ is an interior point of its validity region.

Proof.

Clearly $\overline{\theta}<1$ and hence by Theorem 8.1 ${\cal{L}}(\overline{\theta})$ is an interior point of the parametric-feasible region. We prove by induction that for iteration $\ell$ the base sequence $({\cal{K}},{\cal{J}})^{(\ell)}$ has a validity region with non-empty interior, i.e. there is a $\theta<\overline{\theta}$ big enough, such that the base sequence $({\cal{K}},{\cal{J}})^{(\ell)}$ is optimal for this $\theta$ and ${\cal{L}}(\theta)$ is an interior point of its validity region. Clearly this assumption holds for $\ell=0$. Assume the induction hypothesis for $\ell$ which is optimal for ${\cal{V}}$. Then since validity regions are closed (Theorem S3.7), there exists another validity region with non-empty interior ${\cal{W}}$, such that ${\cal{L}}(\overline{\theta})\in{\cal{W}}$. Furthermore, by the Simplifying Assumption 2 the collision from ${\cal{W}}$ is of type (a)–(f) and hence by Theorem 7.2 base sequence $({\cal{K}},{\cal{J}})^{(\ell+1)}$ that is optimal for ${\cal{W}}$ can be obtained from $({\cal{K}},{\cal{J}})^{(\ell)}$ by a single unique M-CLP pivot. Hence, there exists $\theta>\overline{\theta}$ small enough, such that base sequence $({\cal{K}},{\cal{J}})^{(\ell+1)}$ is optimal for this $\theta$ and ${\cal{L}}(\theta)$ is an interior point of its validity region ${\cal{W}}$. This completes the proof.   ∎

Proposition 8.9**.**

The total number of iterations is bounded by $\left(\!\begin{smallmatrix}4(K+J)\\ 2(K+J)\end{smallmatrix}\!\right)2^{\binom{K+J}{K}}$.

Proof.

Follows from Theorem 8.1.     ∎

9 Illustrative Example Revisited

We now return to the example problem (6) solved in Section 2, and examine it in view of the algorithm of Section 8.

We first list all the admissible bases of the matrix $[A\;I]$. Each of the primal and dual bases is described by the rates of the primal and dual state variables, which determine both the primal and the compatible dual bases. These in fact are the sets ${\cal{K}},{\cal{J}}$.

[TABLE]

For the rate parameters $b,c$ in (18), the objective values of the Rates-LP for these basic solutions are:

[TABLE]

The resulting control rates and state rates for these basic solutions are given in Figure 16.

Note: It is never necessary to calculate all the basic solution of the Rates-LP with matrix $[A\;I]$. The bases which are needed are obtained as necessary when internal pivots are performed, by pivoting from adjacent bases. We list them here to show the rates, and we can do it in this example because it is of such small dimension.

The parametric line of boundary parameters for solving the problem is specified in (19).

We now describe what happens at each of the iterations. We will discuss mainly the boundary pivots. The internal pivots are according to the SCLP algorithm of [36].

Initialization

The initial solution uses zero controls, $U(t)=P(t)=0$, and has states $x(t)=\beta+bt$, $q(t)=-\gamma-ct$. The base sequence includes boundary values and the single internal basis $B_{6}$. We write it as (here we list the internal bases, rather than the $\{{\cal{K}}_{n},{\cal{J}}_{n}\}_{n=1}^{N}$):

[TABLE]

Note that there is a jump in $x(t)$ at $T$: ${\mathbf{x}}^{N}-x^{N}=\lambda-A(U(T)-U(T-))$, and a similar jump in $q(t)$: ${\mathbf{q}}^{0}-q^{0}=-\mu+A^{{\mbox{\tiny\bf\sf T}}}(P(T)-P(T-))$. These jumps are shown in Figure 2. All the state variables in this solution are free.

First iteration

At $\theta^{1}=1/11$ there is a collision of type (e), as $q_{1}(T)={\mathbf{q}}^{0}_{1}$ shrinks to 0. This requires a boundary pivot. No pre-boundary internal pivot is necessary. At this point the boundary values are:

[TABLE]

We can now formulate the Boundary-LP problems (39), (40), and construct the dictionary according to the rules in Section 4.1. In the dictionary, ${\mathbf{q}}^{0}_{1}$ is 0 and needs to leave the dual basis, i.e. $\mathbf{u}^{0}_{1}$ will enter the primal basis. Examination of the dictionary shows that $x_{2}^{\bullet}$ should leave the primal basis, i.e. $\mathbf{p}^{0}_{2}$ will enter the dual basis. Since the test ration is $\neq 0$, this is a boundary pivot of type II.

The boundary dictionaries, before and after the pivot are:

[TABLE]

Following this pivot, in the interval of $(\theta^{(1)},\theta^{(2)})=(1/11,1/6)$ there are two impulse controls: $\mathbf{u}^{0}_{1}$, and $\mathbf{p}^{0}_{2}$, which replaced the boundary values ${\mathbf{q}}^{0}_{1}$ and $x^{0}_{2}$.

2nd to 5th iterations

These are all internal pivots.

At $\theta=1/6$ there is a collision of type (a) at $t_{0}$, where $q^{0}_{1}=q_{1}(T-)$ shrinks to zero. In the internal pivot the internal base sequence is augmented, from $B_{6}$ to $B_{2},B_{6}$. In the collision from the other side, $q_{1}(t)$ becomes free (it is free for $\theta<1/6$ and tied for $\theta>1/6$), so it is a type (d) collision form the other side.

At $\theta=2/9$, there is a collision of type (a) at $t_{1}$, where $x^{1}_{1}=x_{1}(t_{1})$ shrinks to zero. In the internal pivot the internal base sequence is augmented, from $B_{2},B_{6}$ to $B_{2},B_{3},B_{6}$. The collision from the other side is again type (d).

At $\theta=4/9$, there is a collision of type (a) at $t_{1}$, where $q^{1}_{2}=q_{2}(T-t_{1})$ shrinks to zero. In the internal pivot the internal base sequence is augmented, from $B_{2},B_{3},B_{6}$ to $B_{2},B_{1},B_{3},B_{6}$. The collision from the other side is type (d).

At $\theta=5/9$, the interval $(t_{2},t_{3})$ with the rates basis $B_{3}$ shrinks to 0, and $|B_{1}\backslash B_{6}|=|\{u_{1},u_{2}\}\backslash\{\dot{x}_{1},\dot{x}_{2}\}|=2$, so this is a type (b) internal collision. In the internal pivot, $B_{3}$ is deleted from $B_{2},B_{1},B_{3},B_{6}$, and then $B_{4}$ is inserted in its place, so the new internal base sequence is $B_{2},B_{1},B_{4},B_{6}$. The collision form the other side is also type (b).

6th iteration

In the interval $\theta>5/9$ the optimal base sequence is

[TABLE]

All three boundary states in ${\cal{K}}_{0},{\cal{J}}_{N+1}$, namely $x_{1}(t),q_{1}(t),q_{2}(t)$, are now tied.

At $\theta=13/19$, ${\mathbf{x}}^{N}_{2}$ shrinks to zero, which is a collision of type (e), and a boundary pivot is needed. There is no pre-boundary internal pivot. The boundary values at this point are:

[TABLE]

The dictionary for the pivot is:

[TABLE]

The variable that shrunk to 0 is ${\mathbf{x}}^{N}_{2}$, and it will leave the primal basis. We now need to decide on the type of pivot. In this case, the dual boundary values include $q^{\bullet}_{1}=q^{\bullet}_{2}=0$, and both have non-zero pivot elements, so this is a boundary pivot of type I. There is no need to perform a pivot operation on ${\cal{D}}^{(5)}$.

The variable that shrinks to zero on the other side, i.e. from $\theta>13/19$ is the dual variable to ${\mathbf{x}}^{N}_{2}$, which is $\mathbf{p}^{N}_{2}$. The new optimal base sequence is:

[TABLE]

Last iteration

This iteration starts with the optimal base sequence $({\cal{K}},{\cal{J}})^{(6)}$, and this remains optimal for all $13/19<\theta<1$. In the interior of this region there are 4 intervals, and there are three impulse controls, $\mathbf{u}^{0}_{1},\,\mathbf{p}^{0}_{1},\,\mathbf{p}^{N}_{2}$. When $\theta$ reaches 1, a multiple collision occurs, as several items shrink to zero. This is typical, because at that point we have $\lambda=\mu=0$, so the solution of the boundary problem is degenerate. The things that happen at $\theta=1$ are: the interval $(t_{4},T)$ shrinks to zero, the impulse control $\mathbf{p}^{0}_{2}$ shrinks to zero, $x^{N}_{2}$ shrinks to zero, and $q^{N}_{2}$ shrinks to zero.

The final solution of the problem consists of three intervals, and two impulse controls, $\mathbf{u}^{0}_{1}$ at primal time 0, and $\mathbf{p}^{N}_{2}$ at dual time 0. This is the solution displayed in Figure 1.

Appendix A Appendix

A.1 Algorithm under general settings

Under the Simplifying Assumption 2, all collisions were single. We now discuss the algorithm without the simplifying assumption. The main idea is that whenever a multiple collision is discovered in step $\ell$ of the algorithm, going from ${\cal{L}}(\theta_{\ell})$ to ${\cal{L}}(\theta_{\ell+1})$ through an interior of validity region ${\cal{V}}$, one moves to an interior point of ${\cal{V}}$ that is close to the current parametric line, and one restarts from this revised point on a new parametric line leading to ${\cal{L}}(1)$.

The general algorithm works exactly as the simplified algorithm described in Section 8 when collisions at point ${\cal{L}}(\theta_{\ell+1})$ is of type (a)-(f) from both ${\cal{V}}$ and ${\cal{W}}$ sides. However, if collision from the ${\cal{V}}$ or the ${\cal{W}}$ side is a multiple collision there are some additional steps that build a new parametric line ${\cal{L}}^{\prime}(\theta)$.

We extend the collision classification given in Section 6 by the following types:

• •

Pre-$\overline{\theta}$ multiple collision. This is a multiple collision, that is discovered when the solution for the endpoint $\overline{\theta}$ is calculated. It is then seen if the set of values of $H_{\mathbb{P}}^{\ell}$ that shrink to zero

[TABLE]

indicates a multiple collision.

• •

At-$\overline{\theta}$ multiple collision. If there is no pre-$\overline{\theta}$ collision, and there is a collision of types (d)-(f) from the ${\cal{V}}$ side, but in the boundary pivot in which ${\mathbf{v}}$ leaves the boundary basis there are several boundary variables with equal ratio $>0$ that are candidates to enter the basis, this is an at-$\overline{\theta}$ multiple collision. In that case, these candidates to enter are a set of boundary variables that shrink to zero from the ${\cal{W}}$ side as $\theta\searrow\overline{\theta}$. The ratio test, that is a part of boundary pivot discussed in Section 7, performs discovery of this collision.

• •

Post-$\overline{\theta}$ multiple collision. If there is no pre-$\overline{\theta}$ or at-$\overline{\theta}$ multiple collision, then there could be a multiple collision from the ${\cal{W}}$ side in which several local minima of the state variables shrink to zero. This can happens when there is a pivot of type II, and ${\mathbf{w}}^{*}=x_{k}^{\bullet}$ or ${\mathbf{w}}^{*}=q_{j}^{\bullet}$, and it is discovered by the test for multiple collision from ${\cal{W}}$ side discussed in Section 7.

• •

Single collision. This is a collision that occurs under the simplifying assumption, i.e. it is a collision of one of the types (a)-(f) from both the ${\cal{V}}$ and the ${\cal{W}}$ side.

It is easy to see that these collisions cover all possible collisions that may occur along the parametric line.

For all these types of collisions one can choose an interior point ${\cal{L}}(\theta)\in{\cal{V}}$, a direction that is orthogonal to the current parametric line, denoted by $\cal{{\cal{L}}}^{\bot}$, and a step size $\epsilon$ and build a new parametric line ${\cal{L}}^{\prime}(\theta)$ that goes through the two points ${\cal{L}}(\theta)+\epsilon\cal{{\cal{L}}}^{\bot}$ and ${\cal{L}}(1)$ and this line can be chosen to satisfy the following properties:

• •

${\cal{L}}(\theta)+\epsilon\cal{{\cal{L}}}^{\bot}$ is an interior point of ${\cal{V}}$, and the line reaches the boundary of ${\cal{V}}$ at a $\overline{\theta}^{\prime}$.

• •

Parametric line ${\cal{L}}^{\prime}(\theta)$ goes through all the validity regions that were crossed by the parametric line ${\cal{L}}(\theta)$, from the point ${\cal{L}}(0)$ up to the point ${\cal{L}}(\overline{\theta})$.

• •

If there is a pre-$\overline{\theta}$ multiple collision at ${\cal{L}}(\overline{\theta})$ from the ${\cal{V}}$ side, then there is either at-$\overline{\theta}$, post-$\overline{\theta}$ or a single collision at ${\cal{L}}^{\prime}(\overline{\theta}^{\prime})$.

• •

If there is an at-$\overline{\theta}$ multiple collision at ${\cal{L}}(\overline{\theta})$ from the ${\cal{V}}$ side, then there is either post-$\overline{\theta}$ or a single collision at ${\cal{L}}^{\prime}(\overline{\theta}^{\prime})$.

• •

If there is a post-$\overline{\theta}$ multiple collision at ${\cal{L}}(\overline{\theta})$ from the ${\cal{V}}$ side, then there is a single collision at ${\cal{L}}^{\prime}(\overline{\theta}^{\prime})$.

One can see that under such a policy the algorithm cannot cycle and it will still find an optimal solution for the goal set of boundary parameters at the point ${\cal{L}}(1)$ that is shared by all the parametric lines in a finite number of iterations bounded by $\binom{4(K+J)}{2(K+J)}\;2^{\binom{K+J}{K}}$.

We leave the concrete calculations of an interior point ${\cal{L}}^{\prime}(\theta)\in{\cal{V}}$, an orthogonal direction $\cal{{\cal{L}}}^{\bot}$ and a step size $\epsilon$ that satisfy the requirements of this policy out of this paper. We state that this can be done, but it will be more appropriate to deal with the details when one creates an implementation of the algorithm. For further discussion of the requirement of perturbations see [29]

A.2 Discussion of $\lambda$ and $\mu$

The motivation for including $\lambda,\,\mu$ in the M-CLP/M-CLP∗ formulation is that otherwise (i.e. when $\lambda=\mu=0$), solutions in the interior of validity regions may not have $H_{\mathbb{P}}>0$, and validity regions may not have disjoint interiors.

Strong duality and structure properties of M-CLP/M-CLP∗ with $\lambda=\mu=0$ remain valid for solutions of $\lambda,\mu\neq 0$, except for right continuity of solutions at $t=0$. It is possible that optimal solutions now have $U(0-)=0<U(0)<U(0+)$ (or $P(0-)=0<P(0)<P(0+)$).

However, if we restrict the formulation to $\lambda\leq 0$, $\mu\geq 0$, this cannot happen as we state in the next theorem.

Theorem A.1**.**

For feasible M-CLP/M-CLP∗ with $\lambda\leq 0,\mu\geq 0$, there exists a pair of optimal solutions $U_{0}(t),P_{0}(t)$ that are right-continuous at [math].

Proof.

Assume first that the Slater type condition of [30] holds. Imitating the steps leading to Theorem D4.7 and D5.5(iii), M-CLP/M-CLP∗ of (3), (3) possess strongly dual optimal solutions $U_{*}(t),P_{*}(t)$, which are continuous piecewise linear on $(0,T)$. What remains to be checked is whether $U_{*}(t),P_{*}(t)$ are right continuous at [math]. Let $\mathbf{u}_{*}^{0}=U_{*}(0)$ and $\mathbf{u}_{*}^{0+}=U_{*}(0+)-U_{*}(0)$. Consider the following solution of M-CLP:

[TABLE]

We check that $U_{0}(t)$ is a feasible solution of M-CLP at $t=0$. By feasibility of $U_{*}(t)$ we have $AU_{0}(0)=AU_{0}(0+)=AU_{*}(0+)\leq\beta$, and therefore $U_{0}(t)$ is a feasible solution of M-CLP. Next we compare the objective values produced by $U_{*}(t)$ and $U_{0}(t)$.

[TABLE]

Hence, the right-continuous $U_{0}(t)$ is an optimal solution of M-CLP. The construction of $P_{0}(t)$ is similar. Note that $\mu\geq 0$ (and similarly $\lambda\leq 0$) is indeed necessary for the proof. To complete the proof we use the detailed structure of the solution described in [31], where the second boundary equations (S11) have been changed to the equations (30), and use Theorem S4.1 to show that the Slater-type condition is not necessary.     ∎

It is easy to check that all the results of [30, 31] hold for the extended formulation with $\lambda\leq 0$ and $\mu\geq 0$. Further discussion of the motivation and examples are included in [29].

A.3 Proof of Theorem 5.4

Proof.

(i) $\implies$ (iii): Let $\rho$ be an interior point of ${\cal{V}}$ and let $M$, ${\mathbf{R}}$ be the corresponding matrix and r.h.s. of (32) for the base sequence $({\cal{K}}_{n},{\cal{J}}_{n})_{n=0}^{N+1}$ and the point $\rho$. Assume contrary to the theorem that the solution of (32) is not unique.

We now consider the solution of (32), to see which components may be not unique. By Theorem D5.5(iii), all $\tau_{n}$ and $x_{k}^{0},\,k\in{\cal{K}}^{=},\,q_{j}^{N},\,j\in{\cal{J}}^{=}$ are unique. All the boundary values listed in Equations (28) are uniquely determined to be 0. Also, all of $x_{k}(t),\,k\in{\cal{K}}^{=}$ and $\,q_{j}(t),\,j\in{\cal{J}}^{=}$ are uniquely determined by the unique values of the $x_{k}^{0},\,k\in{\cal{K}}^{=},\,q_{j}^{N},\,j\in{\cal{J}}^{=},\,\tau_{n}$. Hence, only components of $\mathbf{u}^{0},\mathbf{u}^{N},\mathbf{p}^{0},\mathbf{p}^{N},{\mathbf{x}}^{N},{\mathbf{q}}^{0}$ and $x_{k}^{0},\,k\in{\cal{K}}^{\uparrow\downarrow},\,q_{j}^{N},\,j\in{\cal{J}}^{\uparrow\downarrow}$ and $x_{k}(t),q_{j}(t),\,k\in{\cal{K}}^{\uparrow\downarrow},j\in{\cal{J}}^{\uparrow\downarrow}$ may be non-unique. By (31), if $x_{k}(t),q_{j}(t),\,k\in{\cal{K}}^{\uparrow\downarrow},j\in{\cal{J}}^{\uparrow\downarrow}$ are non-unique then $x_{k}^{0},\,k\in{\cal{K}}^{\uparrow\downarrow},\,q_{j}^{N},\,j\in{\cal{J}}^{\uparrow\downarrow}$ must be non-unique.

We now note that $\mathbf{u}^{0},\mathbf{u}^{N},\mathbf{p}^{0},\mathbf{p}^{N},{\mathbf{x}}^{N},{\mathbf{q}}^{0}$ and $x_{k}^{0},\,k\in{\cal{K}}^{\uparrow\downarrow},\,q_{j}^{N},\,j\in{\cal{J}}^{\uparrow\downarrow}$ are determined by equations (29) and by equations (30) after substitution of (31). But these equations have $\beta,\gamma,\lambda,\mu$ on the r.h.s..

If the solution is non-unique then $M$ must be singular, and so we must have: $M_{r^{\prime}}=\alpha_{r_{1}}M_{r_{1}}+\dots+\alpha_{r_{m}}M_{r_{m}}$, for which also ${\mathbf{R}}_{r^{\prime}}=\alpha_{r_{1}}{\mathbf{R}}_{r_{1}}+\dots+\alpha_{r_{m}}{\mathbf{R}}_{r_{m}}$ for some of the rows of $M$. But by the above argument, all these rows must have some component of $\beta,\gamma,\lambda,\mu$ on the r.h.s..

Consider $\delta\rho$ in which we change some of the values of $\beta,\gamma,\lambda,\mu$ in such a way that $\delta{\mathbf{R}}_{r^{\prime}}\neq\alpha_{r_{1}}\delta{\mathbf{R}}_{r_{1}}+\dots+\alpha_{r_{m}}\delta{\mathbf{R}}_{r_{m}}$. One can see that for any $\Delta>0$ the equations (32) with r.h.s. ${\mathbf{R}}+\Delta\delta{\mathbf{R}}$ have no solution, and therefore $\rho$ cannot be an interior point. This shows that (i) implies that the solution must be unique.

We now consider the components of $H_{\mathbb{P}}$. Look first at $\tau_{n}$ and $x_{k}(t),\,k\in{\cal{K}}^{=},\,q_{j}(t),\,j\in{\cal{J}}^{=}$. Any of the components of these which belong to $H_{\mathbb{P}}$ can be changed by choosing an appropriate change in $\beta,\gamma$ or $T$. By the above discussion, all the remaining components of $H_{\mathbb{P}}$ can be changed by appropriate choice of $\beta,\gamma,\lambda,\mu$. Assume now that for some element $H_{r}$ of $H_{\mathbb{P}}$, the value is $H_{r}=0$. Then we can choose $\delta\rho$ such that for the corresponding $\delta{\mathbf{R}}$ we have in the solution of $\delta{\mathbf{R}}=M\delta H$ that $\delta H_{r}\neq 0$. But in that case the solution of (32) for $\rho+\Delta\delta\rho$ will have $H_{r}<0$ for all $\Delta>0$ or for all $\Delta<0$, in contradiction to the assumption that $\rho$ is an interior point. This proves that (i) implies that all the components of $H_{\mathbb{P}}$ are positive.

(iii) $\implies$ (ii): Consider $\rho\in{\cal{V}}$ and let $M$, ${\mathbf{R}}$ be the corresponding matrix and r.h.s. of (32) for the base sequence $({\cal{K}}_{n},{\cal{J}}_{n})_{n=0}^{N+1}$ and the point $\rho$. Assume that the solution of (32) for $\rho$ is unique and satisfies $H_{\mathbb{P}}>0$. Hence, for this base sequence the boundary values are unique and strictly complementary slack, which implies the uniqueness of $({\cal{K}}_{n},{\cal{J}}_{n})_{n\in\{0,N+1\}}$. Moreover, all interval lengths are $>0$ and all values of $x^{n},q^{n}$ that are members of $H_{\mathbb{P}}$ are $>0$, and hence at each $t_{n},n=1,\ldots,N-1$ there is exactly one state variable that can leave or enter the basis, and so there cannot be another sequence of bases with zero length intervals between any of $B_{1},\ldots,B_{N}$. Furthermore, at time $T$, none of the variables $x_{k}^{N},\,k\in{\cal{K}}_{N}$ is [math], so no new intervals of length 0 are possible after $B_{N}$, and similarly before $B_{1}$. Hence base sequence $({\cal{K}}_{n},{\cal{J}}_{n})_{n=0}^{N+1}$ is unique.

(ii) $\implies$ (i): Assume $\rho\in{\cal{V}}$ for the unique base sequence $({\cal{K}}_{n},{\cal{J}}_{n})_{n=0}^{N+1}$. Assume contrary to the theorem that $\rho$ is a boundary point of the ${\cal{V}}$. Then, as $\rho$ is an interior point of the parametric-feasible region, from the convexity of the parametric-feasible region and closedness of validity regions (Theorem S3.7) it follows that $\rho$ belongs to a boundary of another validity region ${\cal{V}}^{\prime}$. This contradicts the uniqueness of the $({\cal{K}}_{n},{\cal{J}}_{n})_{n=0}^{N+1}$ and hence $\rho$ is an interior point of the ${\cal{V}}$.   ∎

A.4 Proofs for Section 7

Proof of Theorem 7.1.

The rules for the M-CLP pivot define a unique $({\cal{K}}^{\prime}_{n},{\cal{J}}^{\prime}_{n})_{n=0}^{N^{\prime}+1}$, if we can show the following:

1. (i)

There is no multiple collision discovered from the ${\cal{W}}$ side. 2. (ii)

If the boundary dictionary on the ${\cal{V}}$ side is not unique, then the results of the M-CLP pivot do not depend on the choice of boundary dictionary.

That (i) holds is true by the assumption that the collision from ${\cal{W}}$ is of type (a)–(f) i.e. it is not a multiple collision.

We now show (ii). Consider the case that the leaving variable is a primal variable ${\mathbf{v}}={\mathbf{v}}_{i}=0$ (in row $i$ of the dictionary ${\cal{D}}$). Assume first that for all the dual variables ${\mathbf{w}}_{l}^{*}$ such that ${\mathbf{w}}_{l}^{*}=0$, also ${\hat{{\mathbf{A}}}}_{i,l}=0$. If ${\cal{D}}$ is not a unique dictionary then we can pass from ${\cal{D}}$ to any other dictionary ${\cal{D}}^{\prime\prime}$ by a series of pivots in which a primal variable with value 0 is exchanged for a dual variable of value 0. In that case there must be a pair ${\mathbf{v}}_{k}=0$ and ${\mathbf{w}}_{l}^{*}=0$ such that ${\hat{{\mathbf{A}}}}_{k,l}\neq 0$. Otherwise the dictionary is unique. Consider then the pivot from ${\cal{D}}$ to ${\cal{D}}^{\prime}$ obtained by pivoting on ${\hat{{\mathbf{A}}}}_{k,l}\neq 0$. We assumed that ${\hat{{\mathbf{A}}}}_{i,l}=0$, and as a result, the pivot will not change row $i$. It will also not change the values of the dual variables. Hence the ratio test will yield the same ${\mathbf{w}}^{*}_{r}>0$ that will leave the dual basis, with ${\mathbf{w}}_{r}$ entering the primal basis, for all possible dictionaries.

Assume now that for some ${\mathbf{w}}^{*}_{j}=0$ the element ${\hat{{\mathbf{A}}}}_{i,j}\neq 0$. Then by the above argument, for any other possible dictionary ${\cal{D}}^{\prime\prime}$, there will be some ${\mathbf{w}}^{*}_{l}=0$ such that ${\hat{{\mathbf{A}}}}_{i,l}\neq 0$. But in that case it follows that the pivot under all possible dictionaries will be a pivot of type I. But in that case all that will happen is that ${\mathbf{v}}$ will leave the primal basis, and ${\mathbf{v}}^{*}$ will enter the dual basis, independent of the dictionary chosen.   ∎

Proof of Theorem 7.2.

We distinguish two possibilities.

Case (1): If there are only two optimal base sequences at ${\cal{L}}(\overline{\theta})$, then $({\cal{K}}_{n},{\cal{J}}_{n})_{n=0}^{N+1}$ is optimal in ${\cal{V}}$ and $({\cal{K}}^{\prime}_{n},{\cal{J}}^{\prime}_{n})_{n=0}^{N^{\prime}+1}$ is optimal in ${\cal{W}}$.

• -

Case (2): There is another (or several) optimal base sequence, $({\cal{K}}^{\prime\prime}_{n},{\cal{J}}^{\prime\prime}_{n})_{n=0}^{N^{\prime\prime}+1}$ which is optimal at ${\cal{L}}(\overline{\theta})$. In that case we will show that it has a validity region with an empty interior, so again we have that $({\cal{K}}_{n},{\cal{J}}_{n})_{n=0}^{N+1}$ is optimal in ${\cal{V}}$ and $({\cal{K}}^{\prime}_{n},{\cal{J}}^{\prime}_{n})_{n=0}^{N^{\prime}+1}$ is optimal in ${\cal{W}}$.

We now consider all types of M-CLP pivots, and check if they fall into case (1) or case (2), and in case (2) we show that all other optimal sequences have empty validity regions.

If the pivot is an Internal-SCLP pivot then all the boundary values are positive, so $({\cal{K}}^{\prime}_{n},$ ${\cal{J}}^{\prime}_{n})_{n\in\{0,N^{\prime}+1\}}=({\cal{K}}_{n},{\cal{J}}_{n})_{n\in\{0,N+1\}}$ is unique. Furthermore, by [36], if the collision is not a multiple collision then $({\cal{K}}_{n},{\cal{J}}_{n})_{n=1}^{N}$ and $({\cal{K}}^{\prime}_{n},{\cal{J}}^{\prime}_{n})_{n=1}^{N^{\prime}}$ are the only two optimal base sequences at ${\cal{L}}(\overline{\theta})$. Hence we are in case (1).

If the pivot is a type I boundary pivot and the collision from ${\cal{V}}$ is type (d), then all the boundary values at ${\cal{L}}(\overline{\theta})$ are positive and $({\cal{K}}^{\prime}_{n},{\cal{J}}^{\prime}_{n})_{n\in\{0,N^{\prime}+1\}}=({\cal{K}}_{n},{\cal{J}}_{n})_{n\in\{0,N+1\}}$ is unique, and there are exactly two internal sequences $({\cal{K}}_{n},{\cal{J}}_{n})_{n=1}^{N}$ and $({\cal{K}}^{\prime}_{n},{\cal{J}}^{\prime}_{n})_{n=1}^{N^{\prime}}$ optimal at ${\cal{L}}(\overline{\theta})$. Hence we are in case (1).

If the pivot is a type I boundary pivot and the collision from ${\cal{V}}$ is type (e) or (f), Then the solution at ${\cal{L}}(\overline{\theta})$ is unique, and away from ${\cal{L}}(\overline{\theta})$ either ${\mathbf{v}}>0$ or ${\mathbf{v}}^{*}>0$ determine two unique boundary bases $({\cal{K}}_{n},{\cal{J}}_{n})_{n\in\{0,N+1\}}$ and $({\cal{K}}^{\prime}_{n},{\cal{J}}^{\prime}_{n})_{n\in\{0,N^{\prime}+1\}}$. Further more for each or those there is a unique $({\cal{K}}_{n},{\cal{J}}_{n})_{n=1}^{N}$ and $({\cal{K}}^{\prime}_{n},{\cal{J}}^{\prime}_{n})_{n=1}^{N^{\prime}}$, so again we are in case (1).

Consider now the case of a type II pivot, in which ${\mathbf{v}}=0$ at the collision and is $>0$ on the ${\cal{V}}$ side, and it leaves the primal basis and its dual ${\mathbf{v}}^{*}>0$ enters the dual basis, and ${\mathbf{w}}^{*}>0$ leaves the dual basis and its dual ${\mathbf{w}}=0$ enters the primal basis. Then in $({\cal{K}}_{n},{\cal{J}}_{n})_{n\in\{0,N+1\}}$ we have ${\mathbf{v}}=0$, ${\mathbf{w}}^{*}>0$, and in $({\cal{K}}^{\prime}_{n},{\cal{J}}^{\prime}_{n})_{n\in\{0,N^{\prime}+1\}}$ we have ${\mathbf{v}}^{*}>0$, ${\mathbf{w}}=0$. Apart from these two solutions of the Boundary-LP/LP∗, the only other solutions will be convex combinations of these two solutions, for which there is another base sequence $({\cal{K}}^{\prime\prime}_{n},{\cal{J}}^{\prime\prime}_{n})_{n\in\{0,N^{\prime\prime}+1\}}$ in which the variables ${\mathbf{v}}^{*},{\mathbf{w}}^{*}$ are positive, and at the point ${\cal{L}}(\overline{\theta})$ they can have values on a whole interval. Hence for the base sequence $({\cal{K}}^{\prime\prime}_{n},{\cal{J}}^{\prime\prime}_{n})_{n=0}^{N^{\prime\prime}+1}$ $H_{\mathbb{P}}$ at the point ${\cal{L}}(\overline{\theta})$ is not unique, and hence the matrix $M$ is singular. This implies by Corollary 5.5 that the interior of the validity region of the sequence $({\cal{K}}^{\prime\prime}_{n},{\cal{J}}^{\prime\prime}_{n})_{n=0}^{N^{\prime\prime}+1}$ is empty.   ∎

A.5 Proof of Theorem 8.1

Proof.

(i) By definition ${\cal{L}}(0),{\cal{L}}(1)\in{\cal{F}}$ and the parametric-feasible region is a closed convex polyhedral cone.

(ii) If ${\cal{L}}(0)$ is an interior point of a validity region then it is also an interior point of the parametric-feasible region. By convexity of ${\cal{F}}$ its interior is convex, hence, if point ${\cal{L}}(1)$ is an interior point of the parametric-feasible region then the whole line ${\cal{L}}(\theta)$ belongs to the interior of the parametric-feasible region, otherwise all points of ${\cal{L}}(\theta)$ except ${\cal{L}}(1)$ are interior points of the parametric-feasible region.

(iii) The parametric-feasible region is the union of validity regions of proper base sequences. Hence ${\cal{L}}(\theta)$ is a union of intervals belonging to these validity regions. Because validity regions are convex, no two intervals can belong to the same proper base sequence. By Proposition 5.3 the number of proper base sequences is bounded by $\binom{4(K+J)}{2(K+J)}\;2^{\binom{K+J}{K}}$. Hence the number of intervals is finite. It is then possible to choose a minimal number of proper base sequences which define a minimal number of intervals of positive length.

Note, if part of the parametric line lies on a face of a validity region, then that interval belongs to more than one validity region. In that case this interval will correspond to more than one proper base sequence.     ∎