function [K, f, viol, loc] = hifoo(P, varargin)
%
% HIFOO, A Matlab package for Fixed Order H-infinity and H2 Control
% Stabilization and Performance Optimization for Multiple Plants
%
% K = HIFOO(P) looks for a static output feedback controller that
% stabilizes the plants P{j},j=1,2,... and locally minimizes the max of
% the H-infinity norms of the closed-loop plants.
% K = HIFOO(P, ORDER) looks for an order ORDER controller doing the same.
% K = HIFOO(P, 't') looks for a static output feedback controller that
% stabilizes the plants P{j},j=1,2,... and locally minimizes the max of
% the H2 norms of the closed-loop plants.
% K = HIFOO(P, ORDER, 't') or K = HIFOO(P, 't', ORDER) looks for an
% order ORDER controller doing the same.
% K = HIFOO(P, 'r') looks for a static output feedback controller that
% stabilizes the plants P{j},j=1,2,... and locally minimizes the max of
% reciprocals of the complex stability radii of the closed-loop plants.
% K = HIFOO(P, ORDER, 'r') or K = HIFOO(P, 'r', ORDER) looks for an
% order ORDER controller doing the same.
% K = HIFOO(P, 's') looks for a static output feedback controller that
% locally minimizes the max of the spectral abscissas of the closed-loop
% plants (max of the real parts of their eigenvalues).
% K = HIFOO(P, ORDER, 's') or K = HIFOO(P, 's', ORDER) looks for an
% order ORDER controller doing the same.
% [K, F, VIOL, LOC] = HIFOO(P, ORDER, INIT, FUN, UPPERBND, OPTIONS)
% looks for an order ORDER controller for P{j},j=1,2,... that
% locally solves the following optimization problem:
% min F(K) = max{g_j(K): UPPERBND(j) = +inf}
% subject to the constraints
% g_j(K) <= UPPERBND(j), j=1,2,...,
% where each g_j is one of the supported closed-loop functions
% encoded in FUN, using an initial guess INIT, with options
% given in OPTIONS, returning the controller K, objective value F,
% local optimality certificate LOC and constraint violations VIOL.
%
% Input Parameters:
% P is a cell array of plants, each of which has one of 4 formats
% (if there is only one plant, either P or P{1} may be set to one
% of the first 3 formats):
% - structure with fields
% A,B1,B2,C1,C2,D11,D12,D21,D22
% B can be used instead of B2, C can be used instead of C2,
% D can be used instead of D22.
% If D or D22 is not provided, it is assumed to be zero.
% (B1,C1,D11,D12,D21 are required only if the supported function
% corresponding to the plant is H-infinity performance, and
% are ignored otherwise, OR
% - SS object (as defined in MATLAB's Control System Toolbox),
% encoding A,[B1 B2],[C1; C2],[D11 D12; D21 D22] in standard format
% (for computing H-infinity performance, the index partitioning
% of [B1 B2] and [C1 C2] should be specified in InputGroup.U1,
% InputGroup.U2, OutputGroup.Y1 and OutputGroup.Y2, which are
% set and viewed by "set" and "get" respectively), OR
% - string giving the COMPleib name of the plant
% (see www.compleib.de), OR
% - the single character 'K', equivalent to providing a plant for
% which the closed loop plant equals the controller K.
% ORDER, INIT, FUN, UPPERBND and OPTIONS are optional and may be
% specified in any permutation order
% ORDER: the order of the controller, a nonnegative integer
% (default: 0 (static output feedback))
% INIT: initial guess for controller, has one of 2 formats:
% - structure with fields: a, b, c, d, OR
% - SS object
% When OPTIONS.struct is not used, if the order of the initial guess is
% less than the desired order, the initial guess is augmented to have the
% desired order without increasing the objective value. Thus this routine
% can be called repeatedly to get successively better controllers as the
% order is increased. If the order of the initial guess is greater than
% the desired order, the initial guess is truncated to have the desired
% order but this may increase the objective value.
% If OPTIONS.struct is used, the order of the initial controller must
% equal ORDER.
% (default: INIT is generated randomly)
% FUN: a character or a string, specifying the functions g_j.
% If FUN is a single character, it specifies one supported function
% from the following list that applies to all plants.
% If FUN is a string, it must have length equal to the number of
% plants and it must specify a supported function from the following
% list for each plant; however, plants with UPPERBND(j) = +inf must
% all be associated with the same supported function.
% The supported functions are:
% 'h': +inf if closed loop plant is unstable, otherwise
% H-infinity norm of transfer function from performance input
% to performance output for the closed loop plant
% 't': +inf if closed loop plant is unstable, otherwise H2 norm of
% transfer function from performance input to performance output
% for the closed loop plant
% 'r': +inf if closed loop plant is unstable, otherwise
% reciprocal of complex stability radius of closed loop plant
% 's': spectral abscissa (max(real(eigenvalues))) of closed loop plant
% 'p': pseudospectral abscissa of closed loop plant
% (a value for epsilon must be provided in OPTIONS.epsilon)
% '+': similar to 's', except that spectral abscissa is not optimized:
% hifoo only attempts to find a controller that simultaneously
% stabilizes *all* closed loop plants. Consider using 's' instead.
% (default: 'h')
% UPPERBND: a vector specifying upper bounds on the g_j.
% All plants for which UPPERBND(j) = +inf are unconstrained and
% enter the minimization objective instead.
% Note that if FUN(j) is 'h', 't' or 'r', a stability constraint on
% the closed loop plant for P{j} is implicit as the value of g_j is
% +inf if the closed loop plant is unstable.
% Note: it is allowable to constrain all plants, that is set all
% components of UPPERBND to be finite, unless there is only one plant.
% (default: all +inf)
% OPTIONS: structure with all fields optional:
% OPTIONS.cpumax: quit when cpu time in seconds exceeds this
% (must be a positive number; default: inf)
% OPTIONS.fast:
% 1 to use a fast optimization method (BFGS) only
% 0 to finish optimization with slower method (gradient sampling,
% which may give a better answer and a better output LOC)
% (default: 0)
% OPTIONS.prtlevel: one of 0 (no printing), 1 (minimal, default),
% 2 (includes output from HANSO), 3 (verbose)
% OPTIONS.struct: a structure of matrices with fields
% OPTIONS.struct.a, OPTIONS.struct.b, OPTIONS.struct.c,
% OPTIONS.struct.d consisting of zeros and ones that indicate
% which entries of the controller matrices K.a, K.b, K.c, K.d are
% allowed to be nonzero. For example, if OPTIONS.struct.a(1,1)=0.
% then the controller is restricted to have K.a(1,1)=0, but
% if OPTIONS.struct.a(1,1)=1, K.a(1,1) is free to vary.
% (default: all ones)
% OPTIONS.nrand: number of starting points to use in addition to INIT
% (random perturbations of increasing sizes if INIT is provided)
% (must be a positive integer; default: 3)
% OPTIONS.normtol: termination tolerance (see output parameter LOC)
% (default: 1e-3)
% OPTIONS.evaldist: evaluation distance (see output parameter LOC)
% (default: 1e-3)
% OPTIONS.weight: a vector of positive weights, one for each plant.
% If w = OPTIONS.weight is provided, the problem solved is
% min F(K) = max{w(j)g_j(K): UPPERBND(j) = +inf}
% subject to the constraints
% w(j)*g_j(K) <= UPPERBND(j), j=1,2,...,
% (default: all ones)
% OPTIONS.weightNormK: weight for adding a penalty on the size of
% the controller to the objective function, specifically
% sqrt(||K.a||^2 + ||K.b||^2 + ||K.c||^2 + ||K.d||^2),
% where K.a, ..., K.d are the controller matrices
% (default: 0)
% OPTIONS.augmentHinf: weight for adding reciprocal of complex
% stability radius to H-infinity norm function to avoid closed
% loop plants that are only marginally stabilized: applies to all
% plants for which g_j is the H-infinity norm.
% (default: 0)
% OPTIONS.epsilon: in case that any FUN(j) is 'p', value of epsilon
% defining epsilon-pseudospectral abscissa
% (default: 0.01)
%
% Output parameters
% K: best controller found, has one of 2 formats:
% - structure with 4 fields: a, b, c, d
% - SS object encoding a, b, c, d in SS format
% the format is compatible with the format used for INIT
% (or the format used by P if INIT is not provided)
% F: corresponding value of minimization objective (including any
% terms arising from options.weightNormK and options.augmentHinf,
% but not including any penalty for constraint violations)
% VIOL: a vector specifying constraint violations for the
% plants: VIOL(j) = max(0, g_j(K) - UPPERBND(j))
% (or max(0, w(j)*g_j(K) - UPPERBND(j)) if not unit weights w(j) are
% provided in OPTIONS.weight)
% LOC: local optimality certificate, structure with 2 fields, both set
% to NaN if F is INF or VIOL is nonzero. Otherwise:
% LOC.dnorm: norm of a vector in the convex hull of bundled or
% sampled gradients of the minimization objective,
% evaluated at and near K
% LOC.evaldist: specifies max distance from K at which these
% gradients were evaluated. If this is 0, then LOC.dnorm is
% just the norm of the gradient of the objective at K.
% The smaller LOC.dnorm and LOC.evaldist are, the more likely
% it is that K is an approximate local minimizer.
%
% *Useful Tip*: if the output K is not satisfactory, try calling HIFOO
% again with INIT set to K.
%
% Other software needed
% Required: HANSO 2.0 (Hybrid Algorithm for Non-Smooth Optimization),
% www.cs.nyu.edu/overton/software/hanso/hanso2_0/
% Required when any g_j is 'h' or 'r': either MATLAB's Control System
% Toolbox or SLICOT's linorm (see options.hinfalg, below)
% Required when any g_j is 't': MATLAB's Control System Toolbox
%
% Version 3.0, 2010, see GPL license info below.
% Further documentation is at www.cs.nyu.edu/overton/software/hifoo/
% Marc Millstone, Michael L. Overton, Didier Henrion, Georgia Deaconu,
% Suat Gumussoy, Denis Arzelier
% Send comments/bug reports to Marc Millstone and Michael Overton,
% lastname@cims.nyu.edu, with a subject header containing "hifoo".
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% HIFOO, A Matlab package for Fixed Order H-infinity and H2 control
%% Copyright (C) 2010 Marc Millstone, Michael Overton, Georgia Deaconu
%%
%% This program is free software: you can redistribute it and/or modify
%% it under the terms of the GNU General Public License as published by
%% the Free Software Foundation, either version 3 of the License, or
%% (at your option) any later version.
%%
%% This program is distributed in the hope that it will be useful,
%% but WITHOUT ANY WARRANTY; without even the implied warranty of
%% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
%% GNU General Public License for more details.
%%
%% You should have received a copy of the GNU General Public License
%% along with this program. If not, see .
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%% end of help file%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Additional options which are not intended for typical users, and
% therefore not part of the help info. They may be useful to us or to
% advanced users.
% OPTIONS.rho: a structure used to define the parameter rho
% in the L1 penalty function used to impose the constraints,
% with 3 fields (see below for more information)
% OPTIONS.rho.init: initial value for rho
% OPTIONS.rho.multiply: the amount rho is
% multiplied each time a larger value is tried
% OPTIONS.rho.max: The maximum value allowed for
% rho before quitting
% (default: 100, 100 and 1e10 respectively)
% OPTIONS.maxit: The maximum number of iterations allowed for BFGS
% OPTIONS.hinftol: the tolerance used by the H-infinity norm
% computation (default: 1e-7)
% OPTIONS.hinfalg: a string which selects the algorithm for
% computing the H-infinity norm
% 'matlab' -- Uses the H-infinity norm provided in the Matlab
% Robust Control Toolbox. Experiments show that this
% is slower, but more robust
% 'slicot' -- Uses H-infinity norm, linorm, provided by SLICOT
% (free for non-commercial use). To use this
% functionality, linorm must be in the path.
% Experiments show that this is faster, however,
% a small error may be observed.
% (default ='matlab')
%%%%%%%%%%%%%%%%%%%%%%%%%%%Information about the Method %%%%%%%%%%%%%%%%%%
% There are two phases, Stabilization and Optimization
% STABILIZATION PHASE
% The purpose of this phase is to stablize all plants P{j} for which
% FUN(j) is 'h', 'r' or 't' (H-infinity norm or reciprocal of complex
% stability radius or H2 norm, respectively). Stabilizing these plants,
% we obtain an initial point for which the corresponding g_j is FINITE.
% We apply a combination of BFGS and Gradient Sampling to the max of the
% spectral abscissas (max(real(eig))) of the corresponding closed loop
% plants, starting from INIT and OPTIONS.nrand successively larger random
% perturbations of INIT, or OPTIONS.nrand random starts if INIT is not
% provided, until finding one or more controllers that stabilize these
% plants, or default iteration limits are exceeded.
% If this phase fails, HIFOO terminates with an error message and with
% F or some VIOL(j) (or both) equal to +inf. The stabilization
% phase is skipped if no plants are associated with 'h', 'r' or 't'.
% OPTIMIZATION PHASE
% Using the starting points found by the stabilization phase, or,
% if this was skipped, INIT and OPTIONS.nrand random perturbations,
% we try to optimize the L1 penalty function
% max{g_j(K): UPPERBND(j) = inf} + rho*sum{VIOL(j)}
% where VIOL(j) = max(0, g_j(K) - UPPERBND(j)).
% (If options.weight is not all 1s, weights are included as above;
% if options.weightNormK is not 0, there is an additional term.)
% The initial value for rho is OPTIONS.rho.init. The value of the L1
% penalty function is finite at the starting points and points generated
% with infinite objective value will be rejected by the line search.
% After the minimization is done, a check is made as to whether any
% VIOL(j) is nonzero, and if so, rho is multiplied by OPTIONS.rho.multiply,
% and we try again. This is repeated if necessary until a feasible point
% is found or rho is larger than the maximum value allowed,
% OPTIONS.rho.max, when termination takes place, returning VIOL(j).
% We use only BFGS (not gradient sampling) until we find a feasible point,
% and then call HANSO with the final value for rho to try to
% improve this using gradient sampling (this might possibly require
% increasing the penalty parameter some more.)
% (There is no point in repeating with a larger rho if there are
% no terms in the objective function, that is ALL plants are constrained.)
%
% HANSO 2.0, Hybrid Algorithm for Non-Smooth Optimization, is a
% hybrid of BFGS and gradient sampling. Termination takes place when
% LOC.dnorm <= OPTIONS.normtol with LOC.evaldist approximately <=
% OPTIONS.evaldist, or default iteration limits are exceeded, or
% OPTIONS.cpumax CPU time is exceeded.
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%END OF MAIN COMMENTS%%%%%%%%%%%%%%%%%%%%%%%%