A Generic Solver for Unconstrained Control Problems with Integral Functional Objectives

TL;DR

This paper introduces a versatile and faster solver for unconstrained control problems with integral objectives, significantly improving computational efficiency and demonstrating applicability to classic and extended problems.

Contribution

The paper generalizes and enhances an existing algorithm for unconstrained control problems, achieving faster computation and broader applicability.

Findings

30x faster than previous algorithm on Witsenhausen's counterexample

Successfully applied to additional control problem examples

Potential for extension to constrained control problems

Abstract

We present a generic solver for unconstrained control problems (UCPs) whose objectives take the form of an integral functional of the controllers. The solver generalizes and improves upon the algorithm proposed by Tseng and Tang for the Witsenhausen's counterexample, which provides the best-known results. In essence, we show that minimizing the objective implies minimizing the marginal cost functions almost everywhere, and we perform the latter task pointwisely by the adaptive minimization technique, which speeds up the computation. We implement single-threaded and parallelized versions of the proposed algorithm. Our implementation runs $30 \times$ faster than Tseng and Tang's algorithm on the Witsenhausen's counterexample, and we demonstrate the applicability of the solver and discuss the possible generalization to constrained problems through two more examples.

Equations

$$\min \mathcal{J}[U]$$

$$\mathcal{J}[U]=\int L_m(u_m(y_m),y_m)d{y}_m+{\mathcal{R}}_m[U_{-m}]$$

$$\mathcal{J}[U+ϵ\delta u]=$$

$$+$$

$$\frac{\delta \mathcal{J}[U]}{\delta u_m}(y_m)=0,\frac{\partial }{\partial u_m}\frac{\delta \mathcal{J}[U]}{\delta u_m}(y_m)\ge 0$$

$$u_m(y_m)\leftarrow$$

$$u_m(y_m)=\underset{u\in \mathbb{R}}{\mathrm{argmin}}{L}_m(u,y_m)$$

$$L_m(u_m(y_m),y_m)\ge \Delta +\underset{u\in \mathbb{R}}{\min }{L}_m(u,y_m)$$

$$L_m(u_m(y_m),y_m)=\frac{\partial \mathcal{J}[U]}{\partial y_m}-\frac{\partial {\mathcal{R}}_m[U_{-m}]}{\partial y_m}$$

$$L_m(u,y_m)={\left(\frac{\partial \mathcal{J}[U]}{\partial y_m}-\frac{\partial {\mathcal{R}}_m[U_{-m}]}{\partial y_m}\right)|}_{u_m(y_m)=u}.$$

$$u_m(y_m)=$$

$$=$$

$$C_m(u,y_m)={\frac{\partial \mathcal{J}[U]}{\partial y_m}|}_{u_m(y_m)=u}$$

$$u_m(y_m)=\underset{u\in \mathbb{R}}{\mathrm{argmin}}{C}_m(u,y_m)$$

$$u_m(y_m)\leftarrow \underset{u\in \mathbb{R}}{\mathrm{argmin}}{C}_m(u,y_m).$$

$$\frac{\delta \mathcal{J}[U]}{\delta u_m}(y_m)=$$

$$\frac{\partial }{\partial u_m}\frac{\delta \mathcal{J}[U]}{\delta u_m}(y_m)=$$

$$\left\{u_m(y):y\in B_r(y_m)\right\}$$

$$\left\{u:u\in B_r(u_m(y_m))\right\}.$$

$$u_m(y_m)\leftarrow$$

$$u_m(y_m)\leftarrow y_m\text{for all}{y}_m\in \mathbb{Y}.$$

$$u_m(y_m)\leftarrow$$

$$u_m(y_m)\leftarrow u_m(y_m)-\frac{C_m^{\prime }(u_m(y_m),y_m)}{C_m^{\prime \prime }(u_m(y_m),y_m)}$$

$$C_m^{\prime \prime }(u_m(y_m),y_m)>0$$

$$u_m(y_m)\leftarrow u_m(y_m)-\tau C_m^{\prime }(u_m(y_m),y_m)$$

$$u_m(y_m)\leftarrow$$

$$\min \mathcal{J}[U]=\mathbb{E}\left[k^2{u}_0(y_0{)}^2+x_2^2\right]$$

$$x_1=$$

$$x_2=$$

$$\min \mathcal{J}[U]=\mathbb{E}\left[\lambda u_0(x_0{)}^2+(u_1(x_1)-x_0{)}^2\right]$$

$$y_0=x_0,x_1=u_0(y_0)+w,y_1=x_1,$$

$$\min \mathcal{J}[U]=\mathbb{E}\left[\sum _{m=0}^M\xi u_m(x_m)+\gamma (x_{m+1})\right]$$

Full text

I Introduction

In this work, we focus on the unconstrained control problem (UCP) with the objective

[TABLE]

where $U=\left\{u_{0}(y_{0}),u_{1}(y_{1}),\dots,u_{M-1}(y_{M-1})\right\}$ is the set of controllers $u_{m}:\mathbb{R}\to\mathbb{R}$. We assume that the objective functional $\mathcal{J}[U]$ can be expressed as

[TABLE]

for all $m=0,\cdots,M-1$ (in other words, the functional can be expanded with respect to each $m$), where $L_{m}:\mathbb{R}\times\mathbb{R}\to\mathbb{R}$ and $\mathcal{R}_{m}[U_{-m}]$ is the residual functional depending only on the controllers other than $u_{m}$. UCP is “unconstrained” in the sense that we require the ability to embed all constraints, e.g., system dynamics, in (1). Although omitted, we remark that $L_{m}$ might also depend on $U_{-m}$ and the decomposition of $L_{m}$ and $\mathcal{R}_{m}$ is not necessarily unique.

Finding the optimal controller of a control problem is a daunting task, even when the problem imposes no constraints. The traditional approach to a control problem is to analyze its structure and conclude some useful properties that would help find the optimal controller. However, as demonstrated by the famous Witsenhausen’s counterexample [2], preferred properties, such as linearity, does not hold in general. As a result, deriving the optimal controller becomes craftsmanship relying on keen observations.

Fortunately, we often need a near-optimal controller rather than an exact optimal one in practice. There are two approaches, the analytical and the numerical, to design a near-optimal controller. The former approach examines some specific controller structure based on the problem characteristics. Again, the performance of the design highly depends on the sophisticated understanding of the problem. On the other hand, the numerical approach develops techniques to approximate the problem and obtain good approximations of the optimal controllers. Correspondingly, the main challenges lie on computation efficiency and approximation quality.

It used to be the computation-demanding nature of the numerical methods which impedes the adoption. But luckily, the advancing technologies in the past decades have hugely reshaped the research landscape: Cheaper computation resources and parallelization techniques facilitate the development of image classification [3], machine learning [4, 5], and genomics processing [6]. Increased computation power not only grants us higher efficiency but also enables finer sampling granularity and potentially better approximation quality. Therefore, we argue that numerical methods have great potential in the upcoming computation-rich era.

Although numerical methods could potentially compute faster with plenty of computation resources, its effectiveness plays a central role in getting closer to the optimum. Blindly adopting numerical methods can be ineffective. Accordingly, it is critical to ask how to design numerical methods that are effective for general problems, in particular, for general UCPs.

Our approach to tackling general UCPs evolves from the algorithm proposed for the Witsenhausen’s counterexample in 2017 [1]. The algorithm in [1] outperforms all previous attempts on the well-known Witsenhausen’s counterexample (e.g., [7, 8, 9]), and its mechanism does not depend on the property of the given objective functional. Although the algorithm finds the controllers that result in the record-low cost, it is computationally demanding and requires a special math tool called calculus of variation.

I-A Contribution and Organization

We examine UCP and provide a generic algorithm to find a near-optimal controller numerically. The proposed algorithm generalizes the algorithm in [1] with the following improvements: First, it presents a new angle viewing the local Nash minimizing phase in the algorithm in [1] which does not involve calculus of variation. In summary, our analysis reveals that UCP leads to a per-point marginal cost optimization problem, and the two methods used in the algorithm in [1], local Nash minimizing and local denoising, approach the problem using different candidate sets. We also apply the adaptive minimization technique to speed up the convergence significantly. Furthermore, our proposed algorithm adopts a unified termination criterion applicable to arbitrary UCPs.

We then implement a generic solver based on the proposed algorithm in C++ and demonstrated that it converges $3\times$ faster than the algorithm in [1] on Witsenhausen’s counterexample. Since the proposed algorithm is parallelizable, we enhance the single-threaded C++ implementation using NVIDIA CUDA and observe another $10\times$ computation speed up. We also demonstrate how the solver works on the zero-delay source-channel coding problem, inventory control problem, and $2$-dimensional Witsenhausen’s counterexample. And we open source the tool for future research.

The paper is organized as follows. We first provide a brief overview of the algorithm in [1] in Section II. Section III introduces the ideas of marginal cost functions, local update, partial exhaustion, and adaptive minimization. Those ideas contribute to the design of our solver. We then implement the solver and examine its performance improvement in Section IV. In Section V, we demonstrate that how we can use the solver to approach different problems. Finally, we conclude the paper in Section VI.

I-B Notation

By convention, we denote the state by $x$, control by $u$, observation by $y$, and disturbance by $w$. Let $\mathbb{E}_{A}$ be the expected value with respect to random variable $A$. We omit $A$ when the expectation is taken with respect to all random variables. We denote by $A\sim\mathcal{N}(\mu,\sigma^{2})$ a Gaussian random variable $A$ with mean $\mu$ and variance $\sigma^{2}$ and by $A\sim\mathcal{U}(a,b)$ a uniform random variable distributed over $[a,b]$. Given a number $M$, we slightly abuse the notation to denote $m=0,\dots,M-1$ by $m\in M$. For a multivariate function $F(a,b)$, we introduce the shorthand notation $F^{\prime}(a,b)=\frac{\partial F(a,b)}{\partial a}$ to denote the partial derivative with respect to the first variable.

Given a functional $\mathcal{J}[U]$, we denote by $\frac{\delta\mathcal{J}[U]}{\delta u}(y)$ the functional derivative of $\mathcal{J}[U]$ with respect to the function $u(y)$, which is derived from the Taylor expansion below:

[TABLE]

where $U+\epsilon\delta u$ refers to the set $(U\backslash\left\{u(y)\right\})\cup\left\{u(y)+\epsilon\delta u(y)\right\}$, $\delta u(y)$ is a bounded function for variation, and $O(\epsilon^{3})$ represents the residual terms of order $\epsilon^{3}$ or higher.

II Existing Algorithm and its Limitation

We first explain the algorithm in [1], which attains the state-of-the-art best results for Witsenhausen’s counterexample, in Section II-A. We then discuss the limitations of the algorithm in Section II-B.

II-A Basic Structure

In [1], Witsenhausen’s counterexample is deemed an optimization problem minimizing a given functional. To obtain the best controllers, the algorithm in [1] introduces two main components: local Nash minimizers and local denoising.

II-A1 Local Nash Minimizers

[1] shows that an optimal controller must be a local Nash minimizer. Using calculus of variation, the first order condition (FOC) and the second order condition (SOC)

[TABLE]

are then derived for local Nash minimizers. Combining FOC and SOC, the algorithm in [1] repeats revised Newton’s method to seek for a local Nash minimizer.

II-A2 Local Denoising

The most important observation given by [1] is that finding local Nash minimizers numerically would land in a “noisy” controller. As a result, [1] introduces the idea of “denoising,” i.e., for each $u_{m}$, we denoise for all $y_{m}$ by

[TABLE]

where $B_{r}(a)$ is a ball centering at $a$ with radius $r$.

II-B Limitations

Although the algorithm in [1] is demonstrated effective for Witsenhausen’s counterexample and it is applicable to other similar problems such as inventory control, there are still few issues left by [1]. First, albeit a universal approach to finding local Nash minimizers, calculus of variation is not a simple idea/operation for the people who know little about functional analysis. For local denoising, [1] obtains functions $C_{m}$ by observation. It would be more rigorous to have a standard procedure to derive $C_{m}$. Meanwhile, despite its great performance in terms of the final cost it achieves, the algorithm runs slow in practice. And it is not clear how the termination criterion used in [1] can be easily generalized for arbitrary problems.

III Generic Algorithm Design

In this section, we study UCP and illustrate how to overcome the limitations stated in Section II-B so that we can improve the ideas in [1] to solve UCPs.

III-A Marginal Cost Functions

We start the analysis with Lemma 1, which is a necessary condition for optimal $U$, followed by the derivation of the marginal cost function $C_{m}$.

Lemma 1**.**

If $U$ minimizes $\mathcal{J}$, we have

[TABLE]

almost everywhere for all $m\in M$.

The lemma can be derived from (1) and the derivation is straightforward. Suppose Lemma 1 is not true, there must exist some $\Delta>0$ such that

[TABLE]

over some set $Y_{\Delta}$ with non-zero measure $\mathcal{M}(Y_{\Delta})$. As such, we can set $u_{m}(y_{m})$ as in Lemma 1 for all $y_{m}\in Y_{\Delta}$, which results in the reduction of the functional value by at least $\mathcal{M}(Y_{\Delta})\Delta>0$. However, it leads to a contradiction as $U$ minimizes $\mathcal{J}$.

Lemma 1 relates the optimal $u_{m}(y_{m})$ with $L_{m}(u,y_{m})$. However, a UCP is usually specified by $\mathcal{J}$ and its decomposition to $L_{m}(u,y_{m})$ is in general not unique. We would prefer to base our solver on Lemma 1 with respect to a more deterministic expression than $L_{m}$. Therefore, we derive the marginal cost function $C_{m}$ as follows.

From (1), we know

[TABLE]

which implies that

[TABLE]

Substitute it into Lemma 1, we have

[TABLE]

because $\mathcal{R}_{m}[U_{-m}]$ does not depend on $u$. Thus, we define the marginal cost function $C_{m}$ for each $m\in M$ by

[TABLE]

and we can rephrase Lemma 1 as

Corollary 1**.**

If $U$ minimizes $\mathcal{J}$, we have

[TABLE]

almost everywhere for all $m\in M$.

Now, once we have the objective $\mathcal{J}$ specified, $C_{m}$ can be derived using simple partial derivation. Actually, the definition of $C_{m}$ also matches with the observation in [1].

III-B Local Update and Partial Exhaustion

Corollary 1 is not only a necessary constraint but also a way to improve a non-optimal controller. Essentially, we can always improve a non-optimal solution by

[TABLE]

To do so, we need to find the minimizer of $C_{m}(u,y_{m})$.

To find the minimizer, the most effective method is to exhaust all possible $u\in\mathbb{R}$ and choose the best one. A full exhaustion is rarely practical as the computational cost would be intolerable. Instead, we would prefer computationally economic methods, such as local update and partial exhaustion.

Local update methods include gradient descent, Newton’s method, and their variations. Those methods rely on local information at $y_{m}$, such as derivatives, to update $u_{m}$. Partial exhaustion takes a different approach. It avoids the heavy computational load in a full exhaustion by searching within only a subset of candidate values.

Not surprisingly, these two methods correspond to the two main components in the algorithm in [1]. Since

[TABLE]

finding local Nash minimizers using the revised Newton’s method in [1] is equivalent to minimizing the function $C_{m}(u,y_{m})$ at each $y_{m}$ by a revised Newton’s method.

On the other hand, local denoising is a partial exhaustion method that searches within the candidate set

[TABLE]

to minimize $C_{m}(u,y_{m})$. This is an important feature of the algorithm in [1]: It searches over the domain instead of the range of the function $u_{m}$. Usually, a partial exhaustion method searches the argument value locally, i.e., it considers the candidate set

[TABLE]

Such a partial exhaustion method plays a similar role as a local update method since they both try to update

[TABLE]

for some small $r$. As a result, we do not benefit much from combining the methods together. However, (2) can be a disconnected set when $u_{m}$ is noncontinuous as in Fig. 1, which allows searching and “leaping” to another region.

III-C Adaptive Minimization

Both local update and partial exhaustion methods are adopted alternatively in the algorithm in [1]: In each repetition round, it denoises the controllers after repeating the revised Newton’s method several times. This hybrid strategy converges to the best-known results, but it progresses quite slowly. The reason is that the two methods improve $\mathcal{J}$ in different ways and one may be more effective than the other at different searching phase.

In Fig. 22(a), we consider the Witsenhausen’s counterexample as in [1] and plot the average $\mathcal{J}$ value improvement of the local update method (revised Newton’s method) and the partial exhaustion method (denoising) for a series of repetition rounds, each round consisting of $19$ local update iterations and $1$ partial exhaustion iteration. Fig. 22(a) shows that partial exhaustion performs much better than local update, but local update runs much more times than partial exhaustion. As a result, the overall improvement at each round is low as in Fig. 22(c).

To improve the efficiency, we introduce adaptive minimization. The basic idea is that we will adapt the number of iterations according to the performance of the method. The more effective method gets more iterations to run. Fixing the total iterations to be $20$ per round, we perform adaptive minimization and Fig. 22(b) shows that the average improvements of the two methods are comparable. Furthermore, the overall improvement is boosted by adaptive minimization as shown in Fig. 22(c).

III-D Generic Algorithm

We combine the methods described in Section III-B and Section III-C to construct Algorithm 1 for UCP. In summary, Algorithm 1 approximates the controllers $U$ by step functions over a given sample range. Then, it adaptively minimizes $\mathcal{J}[U]$ using local update and partial exhaustion methods. Adaptive minimization is repeated until some precision factor is met. The whole procedure can be partitioned into five main parts, and we give more detailed descriptions below.

III-D1 Initialization (line 1 – 4)

Given some sampling range $[a,b]\subseteq\mathbb{R}$ and number of samples $d$, we create a linearly spaced vector $\mathbb{Y}\in\mathbb{R}^{d}$ spanning over $[a,b]$ as the domain. Each controller $u_{m}$ is then approximated by a function mapping $\mathbb{Y}$ to $\mathbb{R}$. As suggested in [1], we initialize $u_{m}$ by an identity map:

[TABLE]

Also, we initialize the current objective value $\mathcal{J}_{c}$, improvement $I_{L},I_{P}$, and iterations $N_{L},N_{P}$ for adaptive minimization. We use subscript $L$ to denote the variable for local update and subscript $P$ for partial exhaustion.

III-D2 Main Loop (line 5)

We repeat adaptive minimization in the main loop, and we adopt a unified termination criterion: The loop terminates when the overall improvement in the current round is smaller than the precision $p$. This criterion improves upon the one adopted in [1] as it is valid regardless of the objective $\mathcal{J}$.

III-D3 Local Update (line 6 – 1)

In this phase, we aim to solve

[TABLE]

for all $y_{m}\in\mathbb{Y}$ using some local update methods.

In our implementation, we adopt a modified Newton’s method slightly different from the one in [1]: We only apply Newton’s method

[TABLE]

when

[TABLE]

which implies a local minimum. Otherwise, we apply simple gradient method

[TABLE]

where $\tau$ is the given step size.

III-D4 Partial Exhaustion (line 12 – 1)

We perform partial exhaustion based on the sampled candidate set (2):

[TABLE]

for all $y_{m}\in\mathbb{Y}$. Essentially, it is equivalent to the local denoising procedure with denoising radius $r$. The radius $r$ should be chosen such that it allows the algorithm to examine at least $1$ sample point each direction along $\mathbb{R}$. Meanwhile, $r$ should not be too large as the larger $r$ leads to higher computation load and potentially traps the local update phase around suboptimal local minima. In our implementation, we choose $r$ such that we check $2$ points per direction.

III-D5 Adaptive Minimization (line 18)

We perform a per-round adaptive minimization by updating the iterations according to the improvements in the last round. Starting with fair sharing as in line 1, we calculate the improvements gained from each of the methods by line 1 and 1. In line 18, we allocate the $N$ iterations in the next round proportional to the contribution of each method while ensuring each method will be run at least once.

IV Performance

To demonstrate the performance, we implement Algorithm 1 in two different versions: the single-threaded solver in C++ (denoted by Single) and the parallelized one in NVIDIA CUDA (denoted by Parallel). The solvers are open sourced at [10]. The solvers are run under the number of iterations per round $N=20$ and the precision $p=10^{-10}$.

We compare our solver with the algorithm in [1] on the Witsenhausen’s counterexample, for which the formulation is deferred to Section V-A, under the parameters $k=0.2$ and $\sigma=5$.

For fair comparison, we impose the same termination criterion and total iterations per round for the algorithm in [1]: It performs local Nash minimizing $19$ times followed by $1$ local denoising. All experiments are run on a desktop with Intel Xeon W-2145 (CPU), NVIDIA Quadro P1000 (GPU), and Samsung 970 EVO NVMe M.2 (SSD).

IV-A Parallelization

To reduce the computation time, we examine Algorithm 1 and identify the parts that can be parallelized.

In both the local update and partial exhaustion phases, we revise $u_{m}$ one at a time, which may change $C_{m^{\prime}}$ for some $m^{\prime}\in M,m^{\prime}\neq m$. As a result, such a dependency prevents the computation being done in parallel. On the other hand, updating $u_{m}(y_{m})$ does not affect $C_{m}$, and hence the calculation for each $y_{m}\in\mathbb{Y}$ can be parallelized to reduce the computation time.

IV-B Computation Time

Fig. 3 shows how fast the solvers and the algorithm in [1] converges along the computation time under $14000$ sample points. the algorithm in [1] converges quickly at the beginning since it computes $u_{1}$ directly through a closed form expression. On the contrary, Algorithm 1 does not utilize any closed form expression. It simply searches for $u_{0}$ and $u_{1}$ using local update and partial exhaustion. Although the algorithm in [1] is more effective in computing $u_{1}$, adaptive minimization allows Algorithm 1 to select and use the more effective method. As a result, the single-threaded version converges slower at the beginning, but it improves $\mathcal{J}$ much more effectively than the algorithm in [1] when getting closer to the limit. After parallelization, the convergence speed is further improved by $10\times$, shifting the curve to the left in Fig. 3.

IV-C Scalability

Another property we examine is scalability. Scalability is important for a numerical method as it implies how precise the approximation can be and how large the problem the solver can handle. To demonstrate the scalability, we vary the number of sample points for each controller and measure the computation time needed before reaching the required precision level. In Fig. 4, the single-threaded solver runs about $3\times$ faster than the algorithm in [1]. After parallelization, the computation time is further reduced by $10\times$, and hence it enjoys $30\times$ performance improvement over the algorithm in [1].

V Examples

We apply our solver to three different examples: Witsenhausen’s counterexample (Section V-A), zero-delay source-channel coding (Section V-B), and inventory control (Section V-C). Since inventory control problem imposes additional constraints on the controllers, we discuss how our solver might be generated for the problem with constraints.

V-A Witsenhausen’s Counterexample

Witsenhausen’s counterexample [2] aims to minimize the objective

[TABLE]

subject to the system dynamic

[TABLE]

where $x_{0}\sim\mathcal{N}(0,\sigma^{2})$ and $w\sim\mathcal{N}(0,1)$.

We use our solver to find controllers in Fig. 5. The result is very close to the best-known controller in [1].

V-B Zero-Delay Source-Channel Coding

The zero-delay source-channel coding problem [11] has the objective

[TABLE]

and the system dynamic

[TABLE]

where $x_{0}\sim\mathcal{N}(0,1)$ and $w\sim\mathcal{U}(-1,1)$.

Fig. 6 shows the non-linear controllers found by our solver.

V-C Inventory Control and Constrained Controller

The inventory control problem has the objective

[TABLE]

where $\gamma$ is some cost function attaining the minimum at [math].

The system dynamic is

[TABLE]

where $x_{0}\sim\mathcal{U}(-1,1)$ and $w_{m}\sim\mathcal{U}(-1,1)$ for all $m\in M$.

The inventory control problem imposes an additional constraint: $u_{m}\geq 0$. We can enforce the constraint by performing $u_{m}\leftarrow\max\left\{0,u_{m}\right\}$ after each local update and partial exhaustion. With the same settings as in [1], our solver finds the same controllers as in Fig. 7. It is reported in [1] that this enforcement successfully finds the optimal controllers. Therefore, we conjecture that it is possible to generalize this generic UCP solver for the problems with constraints by projecting the result back to the feasible region after each local update and partial exhaustion.

V-D Multi-Dimensional Controller: $2$-Dimensional Witsenhausen’s Counterexample

We can also apply the solver to multi-dimensional controllers. For example, $2$-dimensional Witsenhausen’s counterexample [12, 13] aims to minimize the objective

[TABLE]

under the same state dynamic as the Witsenhausen’s counterexample [2]. The difference is that the state $x$, controller $u$, and output $y$ are all $2$-dimensional vectors.

We apply our solver with only $160\times 160=25600$ sample points over $y_{0}$ and $y_{1}$ and obtain the results in Fig. 8. The resulting objective value is very close to the best-known one in [13].

VI Conclusion

The paper examines the unconstrained control problems and proposes an effective generic solver which can serve as the benchmark for the future UCP research. On the other hand, a lot of control problems do impose constraints on either the states of the controllers. As a result, it would be of interest to generalize the solver to constrained control problems.