Search

We reformulate the Lanczos tau method for the discretization of time-delay systems in terms of a pencil of operators, allowing for new insights into this approach. As a first main result, we show that, for the choice of a shifted Legendre basis, this method is equivalent to Pad\'e approximation in the frequency domain. We illustrate that Lanczos tau methods straightforwardly give rise to sparse, self nesting discretizations. Equivalence is also demonstrated with pseudospectral collocation, where the non-zero collocation points are chosen as the zeroes of orthogonal polynomials. The importance of such a choice manifests itself in the approximation of the $H^2$-norm, where, under mild conditions, super-geometric convergence is observed and, for a special case, super convergence is proved; both significantly faster than the algebraic convergence reported in previous work., Comment: 20 pages, 4 figures

TDS-CONTROL is an integrated MATLAB package for the analysis and controller-design of linear time-invariant (LTI) dynamical systems with (multiple) discrete delays, supporting both systems of retarded and neutral type. TDS-CONTROL is based on a state-space representations for these TDSs, although functionality is provided to obtain such a formulation from a frequency-domain description of the system. Firstly, the package offers various functionality for analyzing such systems, like methods for computing the spectral abscissa, the H-infinity norm, the pseudospectral abscissa, and the distance to instability. Furthermore, as TDS-CONTROL is designed with neutral time-delay systems in mind, it has the appealing feature that the sensitivity of certain quantities (such as the spectral abscissa) with respect to infinitesimal delay perturbations can explicitly be taken into account. Secondly, TDS-CONTROL also allows to design fixed-order dynamic output feedback controllers. The corresponding controller-design algorithms are based on minimizing the spectral abscissa, the H-infinity norm, or a combination of both with respect to the free controller parameters by solving a non-smooth, non-convex optimization problem. As a strictly negative spectral abscissa is a necessary and sufficient condition for stability, the presented design methods are thus not conservative. It is also possible to impose structure on the controller, enabling the design of decentralized and proportional-integral-derivative (PID) controllers. Furthermore, by allowing the plant to be described in delay descriptor form (i.e., the system's dynamics are given in terms of delay differential algebraic equations), acceleration feedback and Pyragas-type and delay-based controllers can also be considered., Comment: 108 pages, 22 figures, 32 examples

In this work we consider a class of delay eigenvalue problems that admit a spectrum similar to that of a Hamiltonian matrix, in the sense that the spectrum is symmetric with respect to both the real and imaginary axis. More precisely, we present a method to iteratively approximate the eigenvalues of such delay eigenvalue problems closest to a given purely real or imaginary shift, while preserving the symmetries of the spectrum. To this end the presented method exploits the equivalence between the considered delay eigenvalue problem and the eigenvalue problem associated with a linear but infinite-dimensional operator. To compute the eigenvalues closest to the given shift, we apply a specifically chosen shift-invert transformation to this linear operator and compute the eigenvalues with the largest modulus of the new shifted and inverted operator using an (infinite) Arnoldi procedure. The advantage of the chosen shift-invert transformation is that the spectrum of the transformed operator has a "real skew-Hamiltonian"-like structure. Furthermore, it is proven that the Krylov space constructed by applying this operator, satisfies an orthogonality property in terms of a specifically chosen bilinear form. By taking this property into account during the orthogonalization process, it is ensured that even in the presence of rounding errors, the obtained approximation for, e.g., a simple, purely imaginary eigenvalue is simple and purely imaginary. The presented work can thus be seen as an extension of [V. Mehrmann and D. Watkins, "Structure-Preserving Methods for Computing Eigenpairs of Large Sparse Skew-Hamiltonian/Hamiltonian Pencils", SIAM J. Sci. Comput. (22.6), 2001], to the considered class of delay eigenvalue problems. Although the presented method is initially defined on function spaces, it can be implemented using finite dimensional linear algebra operations.

In this paper we solve $m$-parameter eigenvalue problems ($m$EPs), with $m$ any natural number by representing the problem using Tensor-Trains (TT) and designing a method based on this format. $m$EPs typically arise when separation of variables is applied to separable boundary value problems. Often, methods for solving $m$EP are restricted to $m = 3$, due to the fact that, to the best of our knowledge, no available solvers exist for $m>3$ and reasonable size of the involved matrices. In this paper, we prove that computing the eigenvalues of a $m$EP can be recast into computing the eigenvalues of TT-operators. We adapted the algorithm in \cite{Dolgov2014a} for symmetric eigenvalue problems in TT-format to an algorithm for solving generic $m$EPs. This leads to a subspace method whose subspace dimension does not depend on $m$, in contrast to other subspace methods for $m$EPS. This allows us to tackle $m$EPs with $m > 3$ and reasonable size of the matrices. We provide theoretical results and report numerical experiments. The MATLAB code is publicly available., Comment: 17 pages, 14 figures

This manuscript discusses a scalable controller synthesis method for networked systems with a large number of identical subsystems based on the H-infinity control framework. The dynamics of the individual subsystems are described by identical linear time-invariant delay differential equations and the effect of transport and communication delay is explicitly taken into account. The presented method is based on the result that, under a particular assumption on the graph describing the interconnections between the subsystems, the H-infinity norm of the overall system is upper bounded by the robust H-infinity norm of a single subsystem with an additional uncertainty. This work will therefore briefly discuss a recently developed method to compute this last quantity. The resulting controller is then obtained by directly minimizing this upper bound in the controller parameters.

This paper presents the analysis of the stability properties of PID controllers for dynamical systems with multiple state delays, focusing on the mathematical characterization of the potential sensitivity of stability with respect to infinitesimal parametric perturbations. These perturbations originate for instance from neglecting feedback delay, a finite difference approximation of the derivative action, or neglecting fast dynamics. The analysis of these potential sensitivity problems leads us to the introduction of a `robustified' notion of stability called \emph{strong stability}, inspired by the corresponding notion for neutral functional differential equations. We prove that strong stability can be achieved by adding a low-pass filter with a sufficiently large cut-off frequency to the control loop, on the condition that the filter itself does not destabilize the nominal closed-loop system. Throughout the paper, the theoretical results are illustrated by examples that can be analyzed analytically, including, among others, a third-order unstable system where both proportional and derivative control action are necessary for achieving stability, while the regions in the gain parameter-space for stability and strong stability are not identical. Besides the analysis of strong stability, a computational procedure is provided for designing strongly stabilizing PID controllers. Computational case-studies illustrating this design procedure complete the presentation.

In this paper we extend the Residual Arnoldi method for calculating an extreme eigenvalue (e.g. largest real part, dominant,...) to the case where the matrices depend on parameters. The difference between this Arnoldi method and the classical Arnoldi algorithm is that in the former the residual is added to the subspace. We develop a Tensor-Krylov method that applies the Residual Arnoldi method (RA) for a grid of parameter points at the same time. The subspace contains an approximate Krylov space for all these points. Instead of adding the residuals for all parameter values to the subspace we create a low-rank approximation of the matrix consisting of these residuals and add only the column space to the subspace. In order to keep the computations efficient, it is needed to limit the dimension of the subspace and to restart once the subspace has reached the prescribed maximal dimension. The novelty of this approach is twofold. Firstly, we observed that a large error in the low-rank approximations is allowed without slowing down the convergence, which implies that we can do more iterations before restarting. Secondly, we pay particular attention to the way the subspace is restarted, since classical restarting techniques give a too large subspace in our case. We motivate why it is good enough to just keep the approximation of the searched eigenvector. At the end of the paper we extend this algorithm to shift-and-invert Residual Arnoldi method to calculate the eigenvalue close to a shift $\sigma$ for a specific parameter dependency. We provide theoretical results and report numerical experiments. The Matlab code is publicly available., Comment: 20 pages, 5 figures

An eigenvalue based framework is developed for the H-infinity norm analysis and its norm minimization of coupled systems with time-delays, which are naturally described by delay differential algebraic equations (DDAEs). For these equations H-infinity norms are analyzed and their sensitivity with respect to small delay perturbations is studied. Subsequently, numerical methods for the H-infinity norm computation and for designing controllers minimizing the H-infinity norm with a prescribed structure or order, based on a direct optimization approach, are briefly addressed. The effectiveness of the approach is illustrated with a software demo. The chapter concludes by pointing out the similarities with the computation and optimization of characteristic roots of DDAEs., Comment: arXiv admin note: substantial text overlap with arXiv:2003.09252, arXiv:2003.04695, arXiv:2003.05728

H-infinity controllers are frequently used in control theory due to their robust performance and stabilization. Classical H-infinity controller synthesis methods for finite dimensional LTI MIMO plants result in high-order controllers for high-order plants whereas low-order controllers are desired in practice. We design fixed-order H-infinity controllers for a class of time-delay systems based on a non-smooth, non-convex optimization method and a recently developed numerical method for H-infinity norm computations., Comment: Recent Advances in Optimization and its Applications in Engineering, Springer, 2010

We consider the computation of H-infinity norms for Single-Input-Single-Output (SISO) time-delay systems, which are described by delay differential algebraic equations. Unlike the iterative level set methods in the literature, we present a novel numerical method to compute the H-infinity norm. This method requires solving one eigenvalue problem of at most twice the size of the eigenvalue problem in every iteration of a level set method, but in practice often considerably lower. We first show that the computation of extrema of the transfer function can be turned into the computation of the imaginary axis zeros of a transcendental function. We compute these zeros by a predictor-corrector type algorithm. It is known that the H-infinity norm of delay differential algebraic systems, which can model both retarded and neutral type systems, might be sensitive with respect to arbitrarily small delay perturbations. This recently led to the concept of strong H-infinity norms, which explicitly take into account such small delay perturbations. We present a direct numerical method to compute the strong H-infinity norm of SISO time-delay systems. Our algorithm is applicable to the closed-loop system of interconnections (series, parallel, feedback, junctions) of time-delay systems and/or controllers.

We present a numerical method to plot the root locus of Single-Input-Single-Output (SISO) dead-time systems with respect to the controller gain or the system delay. We compute the trajectories of characteristic roots of the closed-loop system on a prescribed complex right half-plane. We calculate the starting, branch and boundary crossing roots of root-locus branches inside the region. We compute the root locus of each characteristic root based on a predictor-corrector type continuation method. To avoid the high sensitivity of roots with respect to the locus parameter in the neighborhood of branch points, the continuation method relies on a natural parameterization of the root-locus trajectory in terms of a distance in the (characteristic root, locus parameter)-space. The method is numerically stable for high order SISO dead-time systems., Comment: arXiv admin note: text overlap with arXiv:2003.04672

We analyze and design H-infinity controllers for general time-delay systems with time-delays in systems' state, inputs and outputs. We allow the designer to choose the order of the controller and to introduce constant time-delays in the controller. The closed-loop system of the plant and the controller is modeled by a system of delay differential algebraic equations (DDAEs). The advantage of the DDAE modeling framework is that any interconnection of systems and controllers prone to various types of delays can be dealt with in a systematic way, without using any elimination technique. We present a predictor-correct algorithm for the H-infinity norm computation of systems described by DDAEs. Instrumental to this we analyze the properties of the H-infinity norm. In particular, we illustrate that it may be sensitive with respect to arbitrarily small delay perturbations. Due to this sensitivity, we introduce the strong H-infinity norm which explicitly takes into account small delay perturbations, inevitable in any practical control application. We present a numerical algorithm to compute the strong H-infinity norm for DDAEs. Using this algorithm and the computation of the gradient of the strong H-infinity norm with respect to the controller parameters, we minimize the strong H-infinity norm of the closed-loop system based on non-smooth, non-convex optimization methods. By this approach, we tune the controller parameters and design H-infinity controllers with a prescribed order or structure., Comment: arXiv admin note: text overlap with arXiv:2003.04695, arXiv:2003.05728

We consider the characterization and computation of H-infinity norms for a class of time-delay systems. It is well known that in the finite dimensional case the H-infinity norm of a transfer function can be computed using the connections between the corresponding singular value curves and the imaginary axis eigenvalues of a Hamiltonian matrix, leading to the established level set methods. We show a similar connection between the transfer function of a time-delay system and the imaginary axis eigenvalues of an infinite dimensional linear operator $\mathcal{L}_{\xi}^N$. Based on this result, we propose a predictor-corrector algorithm for the computation of the H-infinity norm. In the prediction step, a finite-dimensional approximation of the problem, induced by a spectral discretization of the operator $\mathcal{L}_{\xi}$, and an adaptation of the algorithms for finite-dimensional systems, allow to obtain an approximation of the H-infinity norm of the transfer function of the time-delay system. In the next step the approximate results are corrected to the desired accuracy by solving a set of nonlinear equations which are obtained from the reformulation of the eigenvalue problem for the linear infinite-dimensional operator $\mathcal{L}_{\xi}^N$ as a finite dimensional nonlinear eigenvalue problem. These equations can be interpreted as characterizations of peak values in the singular value plot. The effects of the discretization in the predictor step are fully characterized and the choice of the number of discretization points is discussed. The paper concludes with a numerical example and the presentation of the results of extensive benchmarking.

The pseudospectrum of a linear time-invariant system is the set in the complex plane consisting of all the roots of the characteristic equation when the system matrices are subjected to all possible perturbations with a given upper bound. The pseudospectral abscissa is defined as the maximum real part of the characteristic roots in the pseudospectrum and, therefore, it is for instance important from a robust stability point of view. In this paper we present an accurate method for the computation of the pseudospectral abscissa of retarded delay differential equations with discrete pointwise delays. Our approach is based on the connections between the pseudospectrum and the level sets of an appropriately defined complex function. The computation is done in two steps. In the prediction step, an approximation of the pseudospectral abscissa is obtained based on a rational approximation of the characteristic matrix and the application of a bisection algorithm. Each step in this bisection algorithm relies on checking the presence of the imaginary axis eigenvalues of a complex matrix, similar to the delay free case. In the corrector step, the approximate pseudospectral abscissa is corrected to any given accuracy, by solving a set of nonlinear equations that characterize extreme points in the pseudospectrum contours., Comment: arXiv admin note: text overlap with arXiv:2003.04737

An eigenvalue based framework is developed for the stability analysis and stabilization of coupled systems with time-delays, which are naturally described by delay differential algebraic equations. The spectral properties of these equations are analyzed and their stability properties are studied, taking into account the effect of small delay perturbations. Subsequently, numerical methods for stability assessment and for designing stabilizing controllers with a prescribed structure or order, based on a direct optimization approach, are briefly addressed. The effectiveness of the approach is illustrated with a software demo. The paper concludes by pointing out the similarities with the computation and optimization of H-infinity norms., Comment: Delay Systems: Advances in Delays and Dynamics, Springer, 2014

We design fixed-order strong H-infinity controllers for general time-delay systems. The designer chooses the controller order and may introduce constant time-delays in the controller. We represent the closed-loop system of the plant and the controller as delay differential algebraic equations (DDAEs). This representation deals with any interconnection of systems with time-delays without any elimination techniques. We present a numerical algorithm to compute the strong H-infinity norm for DDAEs which is robust to arbitrarily small delay perturbations, unlike the standard H-infinity norm. We optimize the strong H-infinity norm of the closed-loop system based on non-smooth, non-convex optimization methods using this algorithm and the computation of the gradient of the strong H-infinity norm with respect to the controller parameters. We tune the controller parameters and design H-infinity controllers with a prescribed order or structure., Comment: Proceedings of the 18th IFAC World Congress, Milano, Italy, August 28-September 2, 2011. arXiv admin note: text overlap with arXiv:2003.04695

We consider delay differential algebraic equations (DDAEs) to model interconnected systems with time-delays. The DDAE framework does not require any elimination techniques and can directly deal with any interconnection of systems and controllers with time-delays. In this framework, we analyze the properties of the H-infinity norm of systems described by delay differential algebraic equations. We show that the standard H-infinity norm may be sensitive to arbitrarily small delay perturbations. We introduce the strong H-infinity norm which is insensitive to small delay perturbations and describe its properties. We conclude that the strong H-infinity norm is more appropriate in any practical control application compared to the standard H-infinity norm for systems with time-delays whenever there are high-frequency paths in control loops., Comment: Proceedings of the 18th IFAC World Congress, Milano, Italy, August 28-September 2, 2011

We present a numerical method to plot the root-locus of Single-Input-Single-Output (SISO) dead-time systems on a given right half-plane up to a predefined controller gain. We compute the starting and intersection points of root-locus inside the region and we obtain the root-loci of each root based on a predictor-corrector type continuation method. The method is effective for high-order SISO dead-time systems., Comment: Proceedings of the 9th IFAC Workshop on Time-Delay Systems, 2010

The pseudospectra of a linear time-invariant system are the sets in the complex plane consisting of all the roots of the characteristic equation when the system matrices are subjected to all possible perturbations with a given upper bound. The pseudospectral abscissa are defined as the maximum real part of the characteristic roots in the pseudospectra and, therefore, they are for instance important from a robust stability point of view. In this paper we present a numerical method for the computation of the pseudospectral abscissa of retarded delay differential equations with discrete pointwise delays. Our approach is based on the connections between the pseudospectra and the level sets of an appropriately defined complex function. These connections lead us to a bisection algorithm for the computation of the pseudospectral abscissa, where each step relies on checking the presence of imaginary axis eigenvalues of an appropriately defined operator. Because this operator is infinite-dimensional a predictor-corrector approach is taken. In the predictor step the bisection algorithm is applied where the operator is discretized into a matrix, yielding approximations for the pseudospectral abscissa. The effect of the discretization is fully characterized in the paper. In the corrector step, the approximate pseudospectral abscissa are corrected to any given accuracy, by solving a set of nonlinear equations that characterize extreme points in the pseudospectra contours., Comment: arXiv admin note: text overlap with arXiv:2003.03248

In this paper we consider the computation of H-infinity norm of retarded time-delay systems with discrete pointwise state delays. It is well known that in the finite dimensional case H-infinity norm of a system is computed using the connection between the singular values of the transfer function and the imaginary axis eigenvalues of an Hamiltonian matrix. We show a similar connection between the singular values of a transfer function of a time-delay system and the imaginary axis eigenvalues of an infinite dimensional operator $\mathcal{L}_{\xi}$. Using spectral methods, this linear operator is approximated with a matrix. The approximate H-infinity norm of the time-delay system is calculated using the connection between the imaginary eigenvalues of this matrix and the singular values of a finite dimensional approximation of the time-delay system. Finally the approximate results are corrected by solving a set of equations which are obtained from the reformulation of the eigenvalue problem for $\mathcal{L}_{\xi}$ as a finite dimensional nonlinear eigenvalue problem.

This paper examines the robust (strong) H-infinity norm of a linear time-invariant system with discrete delays. The considered system is subject to real-valued, structured, Frobenius norm bounded uncertainties on the coefficient matrices. The robust H-infinity norm is the worst case value of the H-infinity norm over the realisations of the system and hence an important measure of robust performance in control engineering. However this robust H-infinity norm is a fragile measure, as for a particular realization of the uncertainties the H-infinity norm might be sensitive to arbitrarily small perturbations on the delays. Therefore, we introduce the robust strong H-infinity norm, inspired by the notion of strong stability of delay differential equations of neutral type, which takes into account both the perturbations on the system matrices and infinitesimal small delay perturbations. This quantity is a continuous function of the nominal system parameters and delays. The main contribution of this work is the introduction of a relation between this robust strong H-infinity norm and the the pseudo-spectrum of an associated singular delay eigenvalue problem. This relation is subsequently employed in a novel algorithm for computing the robust strong H-infinity norm of uncertain time-delay systems. Both the theoretical results and the algorithm are also generalized to systems with uncertainties on the delays, and systems described by a class of delay differential algebraic equations.

We develop an eigenvalue-based approach for the stability assessment and stabilization of linear systems with multiple delays and periodic coefficient matrices. Delays and period are assumed commensurate numbers, such that the Floquet multipliers can be characterized as eigenvalues of the monodromy operator and by the solutions of a finite-dimensional non-linear eigenvalue problem, where the evaluation of the characteristic matrix involves solving an initial value problem. We demonstrate that such a dual interpretation can be exploited in a two-stage approach for computing dominant Floquet multipliers, where global approximation is combined with local corrections. Correspondingly, we also propose two novel characterizations of left eigenvectors. Finally, from the nonlinear eigenvalue problem formulation, we derive computationally tractable expressions for derivatives of Floquet multipliers with respect to parameters, which are beneficial in the context of stability optimization. Several numerical examples show the efficacy and applicability of the presented results.

In applications of linear algebra including nuclear physics and structural dynamics, there is a need to deal with uncertainty in the matrices. We focus on matrices that depend on a set of parameters $\omega$ and we are interested in the minimal eigenvalue of a matrix pencil $(A,B)$ with $A,B$ symmetric and $B$ positive definite. If $\omega$ can be interpreted as the realisation of random variables, one may be interested in statistical moments of the minimal eigenvalue. In order to obtain statistical moments, we need a fast evaluation of the eigenvalue as a function of $\omega$. Since this is costly for large matrices,we are looking for a small parametrized eigenvalue problem whose minimal eigenvalue makes a small error with the minimal eigenvalue of the large eigenvalue problem. The advantage, in comparison with a global polynomial approximation (on which, e.g., the polynomial chaos approximation relies), is that we do not suffer from the possible non-smoothness of the minimal eigenvalue. The small scale eigenvalue problem is obtained by projection of the large scale problem. Our main contribution is that for constructing the subspace we use multiple eigenvectors as well as derivatives of eigenvectors.We provide theoretical results and document numerical experiments regarding the beneficial effect of adding multiple eigenvectors and derivatives., Comment: 15 pages, 5 figures

A delay Lyapunov matrix corresponding to an exponentially stable system of linear time-invariant delay differential equations can be characterized as the solution of a boundary value problem involving a matrix valued delay differential equation. This boundary value problem can be seen as a natural generalization of the classical Lyapunov matrix equation. We present a general approach for computing delay Lyapunov matrices and H2 norms for systems with multiple discrete delays, whose applicability extends towards problems where the matrices are large and sparse, and the associated positive semidefinite matrix (the ``right-hand side' for the standard Lyapunov equation), has a low rank. In contract to existing methods that are based on solving the boundary value problem directly, our method is grounded in solving standard Lyapunov equations of increased dimensions. It combines several ingredients: i) a spectral discretization of the system of delay equations, ii) a targeted similarity transformation which induces a desired structure and sparsity pattern and, at the same time, favors accurate low rank solutions of the corresponding Lyapunov equation, and iii) a Krylov method for large-scale matrix Lyapunov equations. The structure of the problem is exploited in such a way that the final algorithm does not involve a preliminary discretization step, and provides a fully dynamic construction of approximations of increasing rank. Interpretations in terms of a projection method directly applied to a standard linear infinite-dimensional system equivalent to the original time-delay system are also given. Throughout the paper two didactic examples are used to illustrate the properties of the problem, the challenges and methodological choices, while numerical experiments are presented at the end to illustrate the effectiveness of the algorithm.

This paper is concerned with polynomial approximations of the spectral abscissa function (the supremum of the real parts of the eigenvalues) of a parameterized eigenvalue problem, which are closely related to polynomial chaos approximations if the parameters correspond to realizations of random variables. Unlike in existing works, we highlight the major role of the smoothness properties of the spectral abscissa function. Even if the matrices of the eigenvalue problem are analytic functions of the parameters, the spectral abscissa function may not be everywhere differentiable, even not everywhere Lipschitz continuous, which is related to multiple rightmost eigenvalues or rightmost eigenvalues with multiplicity higher than one. The presented analysis demonstrates that the smoothness properties heavily affect the approximation errors of the Galerkin and collocation-based polynomial approximations, and the numerical errors of the evaluation of coefficients with integration methods. A documentation of the experiments, conducted on the benchmark problems through the software Chebfun, is publicly available., Comment: This is a pre-print of an article published in Numerical Algorithms. The final authenticated version is available online at: https://doi.org/10.1007/s11075-018-00648-9

We consider the minimization or maximization of the $J$th largest eigenvalue of an analytic and Hermitian matrix-valued function, and build on Mengi et al. (2014, SIAM J. Matrix Anal. Appl., 35, 699-724). This work addresses the setting when the matrix-valued function involved is very large. We describe subspace procedures that convert the original problem into a small-scale one by means of orthogonal projections and restrictions to certain subspaces, and that gradually expand these subspaces based on the optimal solutions of small-scale problems. Global convergence and superlinear rate-of-convergence results with respect to the dimensions of the subspaces are presented in the infinite dimensional setting, where the matrix-valued function is replaced by a compact operator depending on parameters. In practice, it suffices to solve eigenvalue optimization problems involving matrices with sizes on the scale of tens, instead of the original problem involving matrices with sizes on the scale of thousands., Comment: 34 pages, 2 figures

We consider a problem in eigenvalue optimization, in particular finding a local minimizer of the spectral abscissa - the value of a parameter that results in the smallest value of the largest real part of the spectrum of a matrix system. This is an important problem for the stabilization of control systems. Many systems require the spectra to lie in the left half plane in order for them to be stable. The optimization problem, however, is difficult to solve because the underlying objective function is nonconvex, nonsmooth, and non-Lipschitz. In addition, local minima tend to correspond to points of non-differentiability and locally non-Lipschitz behavior. We present a sequential linear and quadratic programming algorithm that solves a series of linear or quadratic subproblems formed by linearizing the surfaces corresponding to the largest eigenvalues. We present numerical results comparing the algorithms to the state of the art.