NAG Library Routine Document

e04rpf is a part of the NAG optimization modelling suite and defines bilinear matrix terms either in a new matrix constraint or adds them to an existing linear matrix inequality.

2

Specification

Fortran Interface

Subroutine e04rpf (

handle, nq, qi, qj, dimq, nnzq, nnzqsum, irowq, icolq, q, idblk, ifail)

Integer, Intent (In)	::	nq, qi(nq), qj(nq), dimq, nnzq(nq), nnzqsum, irowq(nnzqsum), icolq(nnzqsum)
Integer, Intent (Inout)	::	idblk, ifail
Real (Kind=nag_wp), Intent (In)	::	q(nnzqsum)
Type (c_ptr), Intent (In)	::	handle

C Header Interface

#include <nagmk26.h>

void	e04rpf_ (void *handle, const Integer nq, const Integer qi[], const Integer qj[], const Integer dimq, const Integer nnzq[], const Integer nnzqsum, const Integer irowq[], const Integer icolq[], const double q[], Integer idblk, Integer ifail)

3

Description

After the initialization routine e04raf has been called, e04rpf may be used to define bilinear matrix terms. It may be used in two ways, either to add to the problem formulation a new bilinear matrix inequality (BMI) which does not have linear terms:

\sum_{i,j = 1}^{n} x_{i} x_{j} Q_{i j} ⪰ 0

(1)

or to extend an existing linear matrix inequality constraint by bilinear terms:

\sum_{i,j = 1}^{n} x_{i} x_{j} Q_{i j}^{k} .

(2)

Here

Q_{i j}^{k}

are

d

d

(sparse) symmetric matrices and

k

, if present, is the number of the existing constraint. This routine will typically be used on semidefinite programming problems with bilinear matrix constraints (BMI-SDP)

\begin{array}{l} \underset{x \in ℝ^{n}}{minimize} & \frac{1}{2} x^{T} H x + c^{T} x & (a) \\ subject to & \sum_{i,j = 1}^{n} x_{i} x_{j} Q_{i j}^{k} + \sum_{i = 1}^{n} x_{i} A_{i}^{k} - A_{0}^{k} ⪰ 0, k = 1, \dots, m_{A} & (b) \\ l_{B} \leq B x \leq u_{B} & (c) \\ l_{x} \leq x \leq u_{x} . & (d) \end{array}

(3)

The routine can be called multiple times to define an additional matrix inequality or to extend an existing one, but it cannot be called twice to extend the same matrix inequality. The argument idblk is used to distinguish whether a new matrix constraint should be added (

idblk = 0

) or if an existing linear matrix inequality should be extended (

idblk > 0

). In the latter case, idblk should be set to

k

, the number of the existing inequality. See e04rnf for details about formulation of linear matrix constraints and their numbering and a further description of idblk. For a generic description of the problem see e04raf. In the further text, the index

k

will be omitted.

3.1

Input data organization

It is expected that only some of the matrices

Q_{i j}

will be nonzero therefore only their index pairs

i, j

are listed in arrays qi and qj. Note that a pair

i, j

should not repeat, i.e., a matrix

Q_{i j}

should not be defined more than once. No particular ordering of pairs

i, j

is expected but other input arrays irowq, icolq, q and nnzq need to respect the chosen order.

Note: the dimension of

Q_{i j}

must respect the size of the linear matrix inequality if they are supposed to expand it (case

idblk > 0

Matrices

Q_{i j}

are symmetric and thus only their upper triangles are passed to the routine. They are stored in sparse coordinate storage format (see Section 2.1.1 in the F11 Chapter Introduction), i.e., every nonzero from the upper triangles is coded as a triplet of row index, column index and the numeric value. All these triplets from all

Q_{i j}

matrices are passed to the routine in three arrays: irowq for row indices, icolq for column indices and q for the values. No particular order of nonzeros within one matrix is enforced but all nonzeros belonging to one

Q_{i j}

matrix need to be stored next to each other. The first

nnzq (1)

nonzeros belong to

Q_{i_{1} j_{1}}

where

i_{1} = qi (1)

j_{1} = qj (1)

, the following

nnzq (2)

nonzeros to the next one given by qi, qj and so on. The array nnzq thus splits arrays irowq, icolq and q into sections so that each section defines one

Q_{i j}

matrix. See Table 1 below. Routines e04rdf and e04rnf use the same data organization so further examples can be found there.

irowq	upper triangle	upper triangle		upper triangle
icolq	nonzeros	nonzeros	$\dots$	nonzeros
q	$\underset{︸}{from Q_{i_{1} j_{1}}}$	$\underset{︸}{from Q_{i_{2} j_{2}}}$		$\underset{︸}{from Q_{i_{nq} j_{nq}}}$
	$nnzq (1)$	$nnzq (2)$		$nnzq (nq)$
	$i_{1} = qi (1)$	$i_{2} = qi (2)$		$i_{nq} = qi (nq)$
	$j_{1} = qj (1)$	$j_{2} = qj (2)$		$j_{nq} = qj (nq)$

Table 1
Coordinate storage format of

Q_{i j}

matrices in input arrays

4

References

Syrmos V L, Abdallah C T, Dorato P and Grigoriadis K (1997) Static output feedback – a survey Journal Automatica (Journal of IFAC) (Volume 33) 2 125–137

5

Arguments

1: $handle$ – Type (c_ptr)Input

On entry: the handle to the problem. It needs to be initialized by e04raf and must not be changed.

2: $nq$ – IntegerInput

On entry: the number of index pairs

i, j

of the nonzero matrices

Q_{i j}

Constraint:

nq > 0

3: $qi (nq)$ – Integer arrayInput

4: $qj (nq)$ – Integer arrayInput

On entry: the index pairs

i, j

of the nonzero matrices

Q_{i j}

in any order.

Constraint:

1 \leq i, j \leq n

where

n

is the number of decision variables in the problem set during the initialization of the handle by e04raf. The pairs do not repeat.

5: $dimq$ – IntegerInput

On entry:

d

, the dimension of matrices

Q_{i j}

Constraints:

$dimq > 0$ ;
if $idblk > 0$ , dimq needs to be identical to the dimension of matrices of the constraint $k$ .

6: $nnzq (nq)$ – Integer arrayInput

On entry: the numbers of nonzero elements in the upper triangles of

Q_{i j}

matrices.

Constraint:

nnzq (i) > 0

7: $nnzqsum$ – IntegerInput

On entry: the dimension of the arrays irowq, icolq and q, at least the total number of all nonzeros in all

Q_{i j}

matrices.

Constraints:

$nnzqsum > 0$ ;
$nnzqsum \geq \sum_{k = 1}^{nq} nnzq (k)$ .

8: $irowq (nnzqsum)$ – Integer arrayInput

9: $icolq (nnzqsum)$ – Integer arrayInput

10: $q (nnzqsum)$ – Real (Kind=nag_wp) arrayInput

On entry: the nonzero elements of the upper triangles of matrices

Q_{i j}

stored in coordinate storage format. The first

nnzq (1)

elements belong to the first

Q_{i_{1} j_{1}}

, the following

nnzq (2)

Q_{i_{2} j_{2}}

, etc.

Constraint:

1 \leq irowq (i) \leq dimq

irowq (i) \leq icolq (i) \leq dimq

11: $idblk$ – IntegerInput/Output

On entry: if

idblk = 0

, a new matrix constraint is created; otherwise

idblk = k > 0

, the number of the existing linear matrix constraint to be expanded with the bilinear terms.

Constraint:

idblk \geq 0

On exit: if

idblk = 0

on entry, the number of the new matrix constraint is returned. By definition, it is the number of the matrix inequalities already defined plus one. Otherwise,

idblk > 0

stays unchanged.

12: $ifail$ – IntegerInput/Output

On entry: ifail must be set to

0

- 1 or 1

. If you are unfamiliar with this argument you should refer to Section 3.4 in How to Use the NAG Library and its Documentation for details.

For environments where it might be inappropriate to halt program execution when an error is detected, the value

- 1 or 1

is recommended. If the output of error messages is undesirable, then the value

1

is recommended. Otherwise, the recommended value is

- 1

. When the value $- 1 or 1$ is used it is essential to test the value of ifail on exit.

On exit:

ifail = 0

unless the routine detects an error or a warning has been flagged (see Section 6).

6

Error Indicators and Warnings

If on entry

ifail = 0

- 1

, explanatory error messages are output on the current error message unit (as defined by x04aaf).

Errors or warnings detected by the routine:

$ifail = 1$: The supplied handle does not define a valid handle to the data structure for the NAG optimization modelling suite. It has not been initialized by e04raf or it has been corrupted.

$ifail = 2$: The problem cannot be modified in this phase any more, the solver has already been called.

$ifail = 3$: On entry, $idblk = 〈value〉$ .
Bilinear terms of the matrix inequality block with the given idblk have already been defined.

$ifail = 4$: On entry, $idblk = 〈value〉$ .
The given idblk does not match with any existing matrix inequality block.
The maximum idblk is currently $〈value〉$ .

On entry, $idblk = 〈value〉$ .
The given idblk refers to a nonexistent matrix inequality block.
No matrix inequalities have been added yet.

$ifail = 5$: On entry, $dimq = 〈value〉$ , $idblk = 〈value〉$ .
The correct dimension of the given idblk is $〈value〉$ .
Constraint: dimq must match the dimension of the block supplied earlier.

$ifail = 6$: On entry, $dimq = 〈value〉$ .
Constraint: $dimq > 0$ .

On entry, $i = 〈value〉$ and $nnzq (i) = 〈value〉$ .
Constraint: $nnzq (i) > 0$ .

On entry, $idblk = 〈value〉$ .
Constraint: $idblk \geq 0$ .

On entry, $nnzqsum = 〈value〉$ and $sum (nnzq) = 〈value〉$ .
Constraint: $nnzqsum \geq sum (nnzq)$ .

On entry, $nq = 〈value〉$ .
Constraint: $nq > 0$ .

$ifail = 8$: On entry, an error occurred in matrix $Q_{i j}$ of index $k = 〈value〉$ , $qi (k) = 〈value〉$ , $qj (k) = 〈value〉$ .
For $j = 〈value〉$ , $icolq (j) = 〈value〉$ and $dimq = 〈value〉$ .
Constraint: $1 \leq icolq (j) \leq dimq$ .

On entry, an error occurred in matrix $Q_{i j}$ of index $k = 〈value〉$ , $qi (k) = 〈value〉$ , $qj (k) = 〈value〉$ .
For $j = 〈value〉$ , $irowq (j) = 〈value〉$ and $dimq = 〈value〉$ .
Constraint: $1 \leq irowq (j) \leq dimq$ .

On entry, an error occurred in matrix $Q_{i j}$ of index $k = 〈value〉$ , $qi (k) = 〈value〉$ , $qj (k) = 〈value〉$ .
For $j = 〈value〉$ , $irowq (j) = 〈value〉$ and $icolq (j) = 〈value〉$ .
Constraint: $irowq (j) \leq icolq (j)$ (elements within the upper triangle).

On entry, an error occurred in matrix $Q_{i j}$ of index $k = 〈value〉$ , $qi (k) = 〈value〉$ , $qj (k) = 〈value〉$ .
More than one element of $Q_{i j}$ has row index $〈value〉$ and column index $〈value〉$ .
Constraint: each element of $Q_{i j}$ must have a unique row and column index.

$ifail = 9$: On entry, index pair with $qi = 〈value〉$ and $qj = 〈value〉$ repeats.
Constraint: each index pair qi, qj must be unique.

On entry, $k = 〈value〉$ , $qi (k) = 〈value〉$ and $n = 〈value〉$ .
Constraint: $1 \leq qi (k) \leq n$ .

On entry, $k = 〈value〉$ , $qj (k) = 〈value〉$ and $n = 〈value〉$ .
Constraint: $1 \leq qj (k) \leq n$ .

$ifail = - 99$: An unexpected error has been triggered by this routine. Please contact NAG.
See Section 3.9 in How to Use the NAG Library and its Documentation for further information.

$ifail = - 399$: Your licence key may have expired or may not have been installed correctly.
See Section 3.8 in How to Use the NAG Library and its Documentation for further information.

$ifail = - 999$: Dynamic memory allocation failed.
See Section 3.7 in How to Use the NAG Library and its Documentation for further information.

7

Accuracy

Not applicable.

8

Parallelism and Performance

e04rpf is not threaded in any implementation.

9

Further Comments

None.

10

Example

This example demonstrates how semidefinite programming can be used in control theory. See also Section 10 in e04raf for links to further examples in the suite.

The problem, from static output feedback (SOF) control Syrmos et al. (1997), solved here is the linear time-invariant (LTI) ‘test’ system

\begin{matrix} \dot{x} = A x + B u \\ y = C x \end{matrix}

(4)

subject to static output feedback

u = K y .

(5)

Here

A \in ℝ^{n \times n}

B \in ℝ^{n \times m}

and

C \in ℝ^{p \times n}

are given matrices,

x \in ℝ^{n}

is the vector of state variables,

u \in ℝ^{m}

is the vector of control inputs,

y \in ℝ^{p}

is the vector of system outputs, and

K \in ℝ^{m \times p}

is the unknown feedback gain matrix.

The problem is to find

K

such that (4) is time-stable when subject to (5), i.e., all eigenvalues of the closed-loop system matrix

A + B K C

belong to the left half-plane. From the Lyapunov stability theory, this holds if and only if there exists a symmetric positive definite matrix

P

such that

{(A + B K C)}^{T} P + P (A + B K C) ≺ 0 .

Hence, by introducing the new variable, the Lyapunov matrix

P

, we can formulate the SOF problem as a feasibility BMI-SDP problem in variables

K

and

P

. As we cannot formulate the problem with sharp matrix inequalities, we can solve the following system instead (note that the objective function is added to bound matrix

P

\begin{array}{l} \underset{K, P}{minimize} & trace (P) \\ subject to & {(A + B K C)}^{T} P + P (A + B K C) ⪯ - I \\ P ⪰ I . \end{array}

(6)

For

n = p = 2

m = 1

A = (\begin{array}{r} - 1 & 2 \\ - 3 & - 4 \end{array}), B = (\begin{matrix} - 1 \\ - 1 \end{matrix}), C = I

and the unknown matrices expressed as

P = (\begin{matrix} x_{1} & x_{2} \\ x_{2} & x_{3} \end{matrix}), K = (\begin{matrix} x_{4} & x_{5} \end{matrix}),

the problem (6) can be rewritten in the form (3) as follows:

\begin{array}{l} \underset{x \in ℝ^{5}}{minimize} & x_{1} + x_{3} \\ subject to & (\begin{array}{r} 2 x_{1} x_{4} + 2 x_{2} x_{4} & x_{1} x_{5} + x_{2} x_{4} + x_{2} x_{5} + x_{3} x_{4} \\ sym. & 2 x_{2} x_{5} + 2 x_{3} x_{5} \end{array}) + \\ (\begin{array}{r} 2 x_{1} + 6 x_{2} & - 2 x_{1} + 5 x_{2} + 3 x_{3} \\ sym. & - 4 x_{2} + 8 x_{3} \end{array}) - I ⪰ 0 \\ (\begin{matrix} x_{1} & x_{2} \\ sym. & x_{3} \end{matrix}) - I ⪰ 0 . \end{array}

This formulation has been stored in a generic BMI-SDP data file which is processed and solved by the example program.

NAG Library Routine Document

e04rpf (handle_set_quadmatineq)

▸▿ Contents

1 Purpose

2 Specification

3 Description

3.1 Input data organization

4 References

5 Arguments

6 Error Indicators and Warnings

7 Accuracy

8 Parallelism and Performance

9 Further Comments

10 Example

10.1 Program Text

10.2 Program Data

10.3 Program Results

1

Purpose

2

Specification

3

Description

3.1

Input data organization

4

References

5

Arguments

6

Error Indicators and Warnings

7

Accuracy

8

Parallelism and Performance

9

Further Comments

10

Example

10.1

Program Text

10.2

Program Data

10.3

Program Results