11.1.4. How to Implement 1-Norm and Infinity-Norm Cost Functions

11.1.4.1. Introduction

In this example we use the system described in the Basic MPC Example, but we will implement non-quadratic costs of the type

\[| \! | Ru_i |\! |_1\]

or

\[| \! | Qx_i |\! |_{\infty}\]

which are sometimes more meaningful for certain applications.

In both cases we will have to introduce slack variables and additional constraints, hence the optimization problem will become more challenging to solve, even if the cost function becomes linear instead of quadratic.

11.1.4.2. 1-norm reformulation

The 1-norm is the absolute sum of a vector, hence a 1-norm penalty on the actuators can be a more meaningful objective when, for instance, the fuel consumption is directly proportional to actuation. The 1-norm also induces sparsity in the solution vector, i.e. a 1-norm cost leads to solutions where actuators are not used at all if possible, which can more accurately represent the objective of minimizing wear in certain applications.

To formulate a 1-norm cost as an optimization problem we introduce one slack variable \(\epsilon_j\) per vector element of \(Ru_i\) (i.e. such that the vector \(\epsilon\) has the same length as the vector \(Ru_i\)) and add it to the polytopic constraints. As a result, the problem

\begin{align*} \text{minimize}& \quad |\! | Ru_i |\! |_1 \\ \text{subject to}& \quad \textrm{constraints} \end{align*}

is transformed into the problem

\begin{align*} \text{minimize}\quad & \sum_{j} \epsilon_j \\ \text{subject to}\quad & \pm Ru_i \leq \epsilon \\ & \textrm{constraints} \end{align*}

The following MATLAB code shows how to model a problem with 1-norm penalties on the actuators and quadratic penalties on the states with FORCESPRO. In particular, note the changes to the cost function and the introduction of polytopic constraints.

%% FORCES multistage form
% assume variable ordering zi = [u{i-1}, x{i}, e{i-1}] for i=1...N

stages = MultistageProblem(N); % get stages struct of length N

for i = 1:N

        % dimension
        stages(i).dims.n = nx+2*nu; % number of stage variables
        stages(i).dims.r = nx;          % number of equality constraints
        stages(i).dims.l = nx+nu;       % number of lower bounds
        stages(i).dims.u = nx+nu;       % number of upper bounds
        stages(i).dims.p = 2*nu;        % number of polytopic constraints

        % cost
        if( i == N )
                stages(i).cost.H = blkdiag(zeros(nu),P,zeros(nu)); % terminal cost (Hessian)
        else
                stages(i).cost.H = blkdiag(zeros(nu),Q,zeros(nu));
        end
        stages(i).cost.f = [zeros(nx+nu,1); ones(nu,1)]; % linear cost terms

        % lower bounds
        stages(i).ineq.b.lbidx = 1:(nu+nx); % lower bound acts on these indices
        stages(i).ineq.b.lb = [umin; xmin]; % lower bound for this stage variable

        % upper bounds
        stages(i).ineq.b.ubidx = 1:(nu+nx); % upper bound acts on these indices
        stages(i).ineq.b.ub = [umax; xmax]; % upper bound for this stage variable

        % polytopic bounds
   stages(i).ineq.p.A = [ R, zeros(nu,nx), -eye(nu); ...
                                                  -R, zeros(nu,nx), -eye(nu)];
        stages(i).ineq.p.b = zeros(2*nu,1);

        % equality constraints
        if( i < N )
                 stages(i).eq.C = [zeros(nx,nu), A, zeros(nx,nu) ];
        end
        if( i>1 )
                 stages(i).eq.c = zeros(nx,1);
        end
        stages(i).eq.D = [B, -eye(nx), zeros(nx,nu)];
end

% RHS of first eq. constr. is a parameter: stages(1).eq.c = -A*x0
params(1) = newParam('minusA_times_x0',1,'eq.c');

You can download the MATLAB code of this example to try it out for yourself by clicking here.

11.1.4.3. \(\infty\)-norm formulation

The \(\infty\)-norm is the maximum absolute value in a vector, hence an \(\infty\)-norm penalty on the states tries to minimize the maximum deviation of any state from the setpoint rather than the combined deviation of all the states in the system.

To formulate an \(\infty\)-norm cost as an optimization problem we need to introduce a single slack variable \(epsilon\) and add polytopic constraints. As a result, the problem

\begin{align*} \text{minimize}& \quad |\! | Qx_i |\! |_{\infty} \\ \text{subject to}& \quad \textrm{constraints} \end{align*}

is transformed into the problem

\begin{align*} \text{minimize}& \quad \epsilon \\ \text{subject to}& \pm Qx_i \leq \mathbf{1}^T \epsilon \\ &\quad \textrm{constraints} \end{align*}

where the vector \(\mathbf{1}=[ 1 \ \ldots \ 1]\) has the same length as the vector \(Qx_i\).

The following MATLAB code shows how to model a problem with ∞-norm penalties on the states and quadratic penalties on the inputs with FORCESPRO. In particular, note the changes to the cost function and the introduction of polytopic constraints. Also note that we only need to add one more variable per stage.

%% FORCES multistage form
% assume variable ordering zi = [u{i-1}, x{i}, e{i-1}] for i=1...N

stages = MultistageProblem(N); % get stages struct of length N

for i = 1:N

        % dimension
        stages(i).dims.n = nx+nu+1; % number of stage variables
        stages(i).dims.r = nx;          % number of equality constraints
        stages(i).dims.l = nx+nu;       % number of lower bounds
        stages(i).dims.u = nx+nu;       % number of upper bounds
        stages(i).dims.p = 2*nx;        % number of polytopic constraints

        % cost
        if( i == N )
                stages(i).cost.H = blkdiag(R,zeros(nx),0); % terminal cost (Hessian)
        else
                stages(i).cost.H = blkdiag(Q,zeros(nx),0);
        end
        stages(i).cost.f = [zeros(nx+nu,1); 1]; % linear cost terms

        % lower bounds
        stages(i).ineq.b.lbidx = 1:(nu+nx); % lower bound acts on these indices
        stages(i).ineq.b.lb = [umin; xmin]; % lower bound for this stage variable

        % upper bounds
        stages(i).ineq.b.ubidx = 1:(nu+nx); % upper bound acts on these indices
        stages(i).ineq.b.ub = [umax; xmax]; % upper bound for this stage variable

        % polytopic bounds
        if( i == N )
           stages(i).ineq.p.A = [ zeros(nx,nu), P, -ones(nx,1); ...
                                                          zeros(nx,nu), -P, -ones(nx,1)];
        else
           stages(i).ineq.p.A = [ zeros(nx,nu), Q, -ones(nx,1); ...
                                                          zeros(nx,nu), -Q, -ones(nx,1)];
        end
        stages(i).ineq.p.b = zeros(2*nx,1);

        % equality constraints
        if( i < N )
                 stages(i).eq.C = [zeros(nx,nu), A, zeros(nx,1)];
        end
        if( i>1 )
                 stages(i).eq.c = zeros(nx,1);
        end
        stages(i).eq.D = [B, -eye(nx), zeros(nx,1)];
end

% RHS of first eq. constr. is a parameter: stages(1).eq.c = -A*x0
params(1) = newParam('minusA_times_x0',1,'eq.c');

You can download the MATLAB code of this example to try it out for yourself by clicking here.