Tutorial¶

First of all, let us import PICOS, CVXOPT, and a pretty printing helper

>>> import picos as pic
>>> import cvxopt as cvx
>>> from pprint import pprint

We now generate some arbitrary data, that we will use in this tutorial.

>>> pairs = [(0,2), (1,4), (1,3), (3,2), (0,4), (2,4)]
>>> b = ([0 ,2 ,0 ,3 ],
...      [1 ,1 ,0 ,5 ],
...      [-1,0 ,2 ,4 ],
...      [0 ,0 ,-2,-1],
...      [1 ,1 ,0 ,0 ])
>>> A = [] # A list of 2×4 matrices.
>>> for i in range(5):
...     A.append(cvx.matrix(range(i-3,i+5),(2,4)))
>>> D = {'Peter': 12,
...      'Bob'  : 4,
...      'Betty': 7,
...      'Elisa': 14}

Let us now create an instance prob of an optimization problem

>>> prob = pic.Problem()

Output settings¶

PICOS makes heavy use of unicode symbols to generate pretty output. If you find that some of these symbols are not available on your terminal, you can call picos.ascii() or picos.latin1() to restrict the charset used.

>>> # Create a dummy symmetric matrix variable to showcase different charsets.
>>> X = pic.Problem().add_variable('X', (2,2), vtype="symmetric")
>>> # Show the default text representation of the SDP constraint X ≽ 0.
>>> print(X >> 0)
X ≽ 0
>>> # Restrict to ISO 8859-1 (aka Latin-1).
>>> pic.latin1()
>>> print(X >> 0)
X » 0
>>> # Restrict to pure ASCII.
>>> pic.ascii()
>>> print(X >> 0)
X >> 0

In addition to complete charset changes, PICOS allows you to fine-tune how certain expressions are formatted. See picos.glyphs for more.

>>> pic.glyphs.psdge.template = "{} - {} is positive semidefinite"
>>> print(X >> 0)
X - 0 is positive semidefinite

For the sake of this tutorial, we return to the full beauty of unicode.

>>> pic.default_charset() # Revert the changes; the same as pic.unicode().

Variables¶

We will now create the variables of our optimization problem. This is done by calling the method add_variable(). This function adds an instance of the class Variable in the dictionary prob.variables, and returns a reference to the freshly added variable. As we will next see, we can use this Variable to form affine and quadratic expressions.

>>> t = prob.add_variable('t',1) #a scalar
>>> x = prob.add_variable('x',4) #a column vector
>>> Y = prob.add_variable('Y',(2,4)) #a matrix
>>> Z = []
>>> for i in range(5):
...     Z.append(prob.add_variable('Z[{0}]'.format(i),(4,2))) # a list of 5 matrices
>>> w = {}
>>> for p in pairs:   #a dictionary of (scalar) binary variables, indexed by our pairs
...     w[p] = prob.add_variable('w[{0}]'.format(p),1 , vtype='binary')

Now, if we try to display a variable, here is what we get:

>>> w[2,4]
<1×1 Binary Variable: w[(2, 4)]>
>>> Y
<2×4 Continuous Variable: Y>

Also note the use of the attributes name, value, size, and vtype:

>>> w[2,4].vtype
'binary'
>>> x.vtype
'continuous'
>>> x.vtype='integer'
>>> x
<4×1 Integer Variable: x>
>>> x.size
(4, 1)
>>> Z[1].value = A[0].T
>>> Z[0].is_valued()
False
>>> Z[1].is_valued()
True
>>> Z[2].name
'Z[2]'

The admissible values for the vtype attribute are documented in add_variable().

Affine Expressions¶

We will now use our variables to create some affine expressions, which are stored as instance of the class AffinExp, and will be the core to define an optimization problem. Most python operators have been overloaded to work with instances of AffinExp (a list of available overloaded operators can be found in the doc of AffinExp). For example, you can form the sum of two variables by writing:

>>> Z[0]+Z[3]
<4×2 Affine Expression: Z[0] + Z[3]>

The transposition of an affine expression is done by appending .T:

>>> x
<4×1 Integer Variable: x>
>>> x.T
<1×4 Affine Expression: xᵀ>

Parameters as constant affine expressions¶

It is also possible to form affine expressions by using parameters stored in data structures such as a list or a cvxopt matrix (In fact, any type that is recognizable by the function retrieve_matrix).

>>> x + b[0]
<4×1 Affine Expression: x + [4×1]>
>>> x.T + b[0]
<1×4 Affine Expression: xᵀ + [1×4]>
>>> A[0] * Z[0] + A[4] * Z[4]
<2×2 Affine Expression: [2×4]·Z[0] + [2×4]·Z[4]>

In the above example, you see that the list b[0] was correctly converted into a $4 \times 1$ vector in the first expression, and into a $1 \times 4$ vector in the second one. This is because the overloaded operators always try to convert the data into matrices of the appropriate size.

If you want to have better-looking string representations of your affine expressions, you will need to convert the parameters into constant affine expressions. This can be done thanks to the function new_param():

>>> A = pic.new_param('A',A)              #this creates a list of constant affine expressions [A[0],...,A[4]]
>>> b = pic.new_param('b',b)              #this creates a list of constant affine expressions [b[0],...,b[4]]
>>> D = pic.new_param('D',D)              #this creates a dictionary of constant AffExpr, indexed by 'Peter', 'Bob', ...
>>> alpha = pic.new_param('alpha',12)     #a scalar parameter

>>> alpha
<1×1 Affine Expression: alpha>
>>> D['Betty']
<1×1 Affine Expression: D[Betty]>
>>> pprint(b)
[<4×1 Affine Expression: b[0]>,
 <4×1 Affine Expression: b[1]>,
 <4×1 Affine Expression: b[2]>,
 <4×1 Affine Expression: b[3]>,
 <4×1 Affine Expression: b[4]>]
>>> print(b[0])
[ 0.00e+00]
[ 2.00e+00]
[ 0.00e+00]
[ 3.00e+00]
<BLANKLINE>

The above example also illustrates that when a valued affine expression is printed, then its value is displayed instead of a symbolic description. Note that the constant affine expressions, as b[0] in the above example, are always valued. To assign a value to a non-constant AffinExp, you must set the value property of every variable involved in the affine expression.

>>> x_minus_1 = x - 1
>>> x_minus_1                           #note that 1 was recognized as the (4x1)-vector with all ones
<4×1 Affine Expression: x - [1]>
>>> print(x_minus_1)
x - [1]
>>> x_minus_1.is_valued()
False
>>> x.value = [0,1,2,-1]
>>> x_minus_1.is_valued()
True
>>> print(x_minus_1)
[-1.00e+00]
[ 0.00e+00]
[ 1.00e+00]
[-2.00e+00]
<BLANKLINE>

We also point out that new_param() converts lists into vectors and lists of lists into matrices (given in row major order). In contrast, tuples are converted into list of affine expressions:

>>> pic.new_param('vect',[1,2,3])                        # [1,2,3] is converted into a vector of dimension 3
<3×1 Affine Expression: vect>
>>> pic.new_param('mat',[[1,2,3],[4,5,6]])               # [[1,2,3],[4,5,6]] is converted into a (2x3)-matrix
<2×3 Affine Expression: mat>
>>> pic.new_param('list_of_scalars',(1,2,3))             # (1,2,3) is converted into a list of 3 scalar parameters #doctest: +NORMALIZE_WHITESPACE
[<1×1 Affine Expression: list_of_scalars[0]>,
 <1×1 Affine Expression: list_of_scalars[1]>,
 <1×1 Affine Expression: list_of_scalars[2]>]
>>> pic.new_param('list_of_vectors',([1,2,3],[4,5,6]))   # ([1,2,3],[4,5,6]) is converted into a list of 2 vector parameters #doctest: +NORMALIZE_WHITESPACE
[<3×1 Affine Expression: list_of_vectors[0]>,
 <3×1 Affine Expression: list_of_vectors[1]>]

Overloaded operators¶

OK, so now we have some variables (t, x, w, Y, and Z) and some parameters (A, b, D and alpha). Let us create some affine expressions with them.

>>> A[0] * Z[0]                              #left multiplication
<2×2 Affine Expression: A[0]·Z[0]>
>>> Z[0] * A[0]                              #right multiplication
<4×4 Affine Expression: Z[0]·A[0]>
>>> A[1] * Z[0] * A[2]                       #left and right multiplication
<2×4 Affine Expression: A[1]·Z[0]·A[2]>
>>> alpha*Y                                  #scalar multiplication
<2×4 Affine Expression: alpha·Y>
>>> t/b[1][3] - D['Bob']                     #division by a scalar and substraction
<1×1 Affine Expression: t/b[1][3] - D[Bob]>
>>> ( b[2] | x )                             #dot product
<1×1 Affine Expression: ⟨b[2], x⟩>
>>> ( A[3] | Y )                             #generalized dot product for matrices: (A|B)=trace(A*B.T)
<1×1 Affine Expression: ⟨A[3], Y⟩>
>>> b[1]^x                                   #hadamard (element-wise) product
<4×1 Affine Expression: b[1]⊙x>

We can also take some subelements of affine expressions, by using the standard syntax of python slices:

>>> b[1][1:3]                                #2d and 3rd elements of b[1]
<2×1 Affine Expression: b[1][1:3]>
>>> Y[1,:]                                   #2d row of Y
<1×4 Affine Expression: Y[1,:]>
>>> x[-1]                                    #last element of x
<1×1 Affine Expression: x[-1]>
>>> A[2][:,1:3]*Y[:,-2::-2]                  #extended slicing with (negative) steps is allowed
<2×2 Affine Expression: A[2][:,1:3]·Y[:,-2::-2]>

In the last example, we keep only the second and third columns of A[2], and the columns of Y with an even index, considered in the reverse order. To concatenate affine expressions, the operators // (vertical concatenation) and & (horizontal concatenation) have been overloaded:

>>> (b[1] & b[2] & x & A[0].T*A[0]*x) // x.T
<5×4 Affine Expression: [b[1], b[2], x, A[0]ᵀ·A[0]·x; xᵀ]>

When a scalar is added/substracted to a matrix or a vector, we interprete it as an elementwise addition of the scalar to every element of the matrix or vector.

>>> 5*x - alpha
<4×1 Affine Expression: 5·x - [alpha]>

Summing Affine Expressions¶

You can take the advantage of python syntax to create sums of affine expressions:

>>> sum([A[i]*Z[i] for i in range(5)])
<2×2 Affine Expression: A[0]·Z[0] + A[1]·Z[1] + A[2]·Z[2] + A[3]·Z[3] + A[4]·Z[4]>

This works, but you might have very long string representations if there are a lot of summands. So you’d better use the function picos.sum()):

>>> pic.sum([A[i]*Z[i] for i in range(5)])
<2×2 Affine Expression: ∑(A[i]·Z[i] : i ∈ [0…4])>

Objective function¶

The objective function of the problem can be defined with the function set_objective(). Its first argument should be 'max', 'min' or 'find' (for feasibility problems), and the second argument should be a scalar expression:

>>> prob.set_objective('max',( A[0] | Y )-t)
>>> print(prob) #doctest: +NORMALIZE_WHITESPACE
---------------------
optimization problem (Mixed-Integer Unconstrained Problem):
59 variables, 0 affine constraints
<BLANKLINE>
Y   : (2, 4), continuous
Z   : list of 5 variables, (4, 2), continuous
t   : (1, 1), continuous
w   : dict of 6 variables, (1, 1), binary
x   : (4, 1), integer
<BLANKLINE>
    maximize ⟨A[0], Y⟩ - t
such that
  []
---------------------

With this example, you see what happens when a problem is printed: the list of optimization variables is displayed, then the objective function and finally a list of constraints (in the case above, there is no constraint).

Norm of an affine Expression¶

The norm of an affine expression is an overload of the abs() function. If x is an affine expression, abs(x) is its Euclidean norm $\sqrt{x^T x}$ .

>>> abs(x)
<Norm of a 4×1 Expression: ‖x‖>

In the case where the affine expression is a matrix, abs() returns its Frobenius norm, defined as $\Vert M \Vert_F := \sqrt{\operatorname{trace} (M^T M)}$ .

>>> abs(Z[1]-2*A[0].T)
<Norm of a 4×2 Expression: ‖Z[1] - 2·A[0]ᵀ‖>

Note that the absolute value of a scalar expression is stored as a norm:

>>> abs(t)
<Norm of a 1×1 Expression: ‖t‖>

However, a scalar constraint of the form $|a^T x + b| \leq c^T x + d$ is handled as two linear constraints by PICOS, and so a problem with the latter constraint can be solved even if you do not have a SOCP solver available.

It is also possible to use other $L_p-$ norms in picos, cf. this paragraph .

Quadratic Expressions¶

Quadratic expressions can be formed in several ways:

>>> t**2 - x[1]*x[2] + 2*t - alpha                       #sum of linear and quadratic terms
<Quadratic Expression: t² - x[1]·x[2] + 2·t - alpha>
>>> (x[1]-2) * (t+4)                                     #product of two affine expressions
<Quadratic Expression: (x[1] - 2)·(t + 4)>
>>> Y[0,:]*x                                             #Row vector multiplied by column vector
<Quadratic Expression: Y[0,:]·x>
>>> (x +2 | Z[1][:,1])                                   #scalar product of affine expressions
<Quadratic Expression: ⟨x + [2], Z[1][:,1]⟩>
>>> abs(x)**2                                            #recall that abs(x) is the euclidean norm of x
<Quadratic Expression: ‖x‖²>
>>> (t & alpha) * A[1] * x                               #quadratic form
<Quadratic Expression: [t, alpha]·A[1]·x>

It is not possible (yet) to make a multidimensional quadratic expression.

Constraints¶

A constraint takes the form of two expressions separated by a relation operator.

Linear (in)equalities¶

Linear (in)equalities are understood elementwise. The strict operators < and > denote weak inequalities (less or equal than and larger or equal than). For example:

>>> (1|x) < 2                                                        #sum of the x[i] less or equal than 2
<1×1 Affine Constraint: ⟨[1], x⟩ ≤ 2>
>>> Z[0] * A[0] > b[1]*b[2].T                                        #A 4x4-elementwise inequality
<4×4 Affine Constraint: Z[0]·A[0] ≥ b[1]·b[2]ᵀ>
>>> pic.sum([A[i]*Z[i] for i in range(5)]) == 0            #A 2x2 equality. The RHS is the all-zero matrix
<2×2 Affine Constraint: ∑(A[i]·Z[i] : i ∈ [0…4]) = 0>

Constraints can be added in the problem with the function add_constraint():

>>> constraints = []
>>> for i in range(1,5):
...     constraints.append(prob.add_constraint(Z[i]==Z[i-1]+Y.T))
>>> print(prob) #doctest: +NORMALIZE_WHITESPACE
---------------------
optimization problem (MILP):
59 variables, 32 affine constraints
<BLANKLINE>
Y   : (2, 4), continuous
Z   : list of 5 variables, (4, 2), continuous
t   : (1, 1), continuous
w   : dict of 6 variables, (1, 1), binary
x   : (4, 1), integer
<BLANKLINE>
    maximize ⟨A[0], Y⟩ - t
such that
  Z[1] = Z[0] + Yᵀ
  Z[2] = Z[1] + Yᵀ
  Z[3] = Z[2] + Yᵀ
  Z[4] = Z[3] + Yᵀ
---------------------

The constraints of the problem can be accessed via the reference returned by add_constraint() or with the function get_constraint():

>>> prob.get_constraint(2)
<4×2 Affine Constraint: Z[3] = Z[2] + Yᵀ>
>>> constraints[2] is prob.get_constraint(2)
True

Groupping constraints¶

In order to have a more compact string representation of the problem, it is advised to use the function add_list_of_constraints(), which works similarly as the function sum().

>>> # We first remove all constraints previously added.
>>> prob.remove_all_constraints()
>>> # Then we add a single constraint.
>>> C = prob.add_constraint(Y > 0)
>>> # And then the same list of constraints as above.
>>> CL = prob.add_list_of_constraints(
...     [Z[i] == Z[i-1] + Y.T for i in range(1,5)])
>>> print(prob) #doctest: +NORMALIZE_WHITESPACE
---------------------
optimization problem (MILP):
59 variables, 40 affine constraints
<BLANKLINE>
Y   : (2, 4), continuous
Z   : list of 5 variables, (4, 2), continuous
t   : (1, 1), continuous
w   : dict of 6 variables, (1, 1), binary
x   : (4, 1), integer
<BLANKLINE>
    maximize ⟨A[0], Y⟩ - t
such that
  Y ≥ 0
  Z[i+1] = Z[i] + Yᵀ ∀ i ∈ [0…3]
---------------------

Now, the constraint Z[3] == Z[2] + Y.T, which has been entered in 4th position, can either be accessed by prob.get_constraint(3) (3 because constraints are numbered from 0), or by

>>> prob.get_constraint((1,2))
<4×2 Affine Constraint: Z[3] = Z[2] + Yᵀ>

where (1,2) means the 3rd constraint of the 2d group of constraints, with zero-based numbering.

Similarly, the constraint Y ≥ 0 can be accessed by prob.get_constraint(0) (first constraint), prob.get_constraint((0,0)) (first constraint of the first group), or prob.get_constraint((0,)) (unique constraint of the first group).

Removing constraints¶

It can be useful to remove some constraints, especially for dynamic approaches such as column generation. Re-creating an instance from scratch after each iteration would be inefficient. Instead, PICOS allows one to modify the solver instance and to re-solve it on the fly (for mosek, cplex and gurobi). It is not possible to change directly a constraint, but you can delete a constraint from model, and then re-add the modified constraint.

The next code shows an example in which a variable x2 is added to the model, which appears as +3*x2 in the objective function and as +x2 in the LHS of a constraint.

>>> prb = pic.Problem()
>>> x1 = prb.add_variable('x1',1)
>>> lhs = 2*x1
>>> obj = 5*x1
>>> cons = prb.add_constraint(lhs <= 1)
>>> sol = prb.maximize(obj,verbose=0,solver='cvxopt')
>>> #--------------------------------------
>>> #at this place, the user can use his favorite method to solve the 'pricing problem'.
>>> #Let us assume that this phase suggests to add a new variable 'x2' in the model,
>>> #which appears as `+3*x2` in the objective function, and as `+x2` in the LHS of the constraint.
>>> x2 = prb.add_variable('x2',1)
>>> lhs += x2
>>> obj += 3*x2
>>> cons.delete()
>>> newcons = prb.add_constraint(lhs <= 1)
>>> print(prb) #doctest: +NORMALIZE_WHITESPACE
---------------------
optimization problem (LP):
2 variables, 1 affine constraints
<BLANKLINE>
x1  : (1, 1), continuous
x2  : (1, 1), continuous
<BLANKLINE>
    maximize 5·x1 + 3·x2
such that
  2·x1 + x2 ≤ 1
---------------------

Flow constraints in Graphs¶

Flow constraints in graphs are entered using a Networkx Graph. The following example finds a (trivial) maximal flow from 'S' to 'T' in G.

>>> import networkx as nx
>>> G = nx.DiGraph()
>>> G.add_edge('S','A', capacity=1)
>>> G.add_edge('A','B', capacity=1)
>>> G.add_edge('B','T', capacity=1)
>>> pb = pic.Problem()
>>> # Adding the flow variables
>>> f={}
>>> for e in G.edges():
...     f[e]=pb.add_variable('f[{0}]'.format(e),1)
>>> # A variable for the value of the flow
>>> F = pb.add_variable('F',1)
>>> # Creating the flow constraint
>>> flowCons = pb.add_constraint(pic.flow_Constraint(
...     G, f, 'S', 'T', F, capacity='capacity', graphName='G'))
>>> pb.set_objective('max',F)
>>> sol = pb.solve(verbose=0,solver='cvxopt')
>>> flow = pic.tools.eval_dict(f)

Picos allows you to define single source - multiple sinks problems. You can use the same syntax as for a single source - single sink problems. Just add a list of sinks and a list of flows instead.

import picos as pic
import networkx as nx

G=nx.DiGraph()
G.add_edge('S','A', capacity=2); G.add_edge('S','B', capacity=2)
G.add_edge('A','T1', capacity=2); G.add_edge('B','T2', capacity=2)

pbMultipleSinks=pic.Problem()
# Flow variable
f={}
for e in G.edges():
    f[e]=pbMultipleSinks.add_variable('f[{0}]'.format(e),1)

# Flow value
F1=pbMultipleSinks.add_variable('F1',1)
F2=pbMultipleSinks.add_variable('F2',1)

flowCons = pic.flow_Constraint(
    G, f, source='S', sink=['T1','T2'], capacity='capacity',
    flow_value=[F1, F2], graphName='G')

pbMultipleSinks.add_constraint(flowCons)
pbMultipleSinks.set_objective('max',F1+F2)

# Solve the problem
pbMultipleSinks.solve(verbose=0,solver='cvxopt')

print(pbMultipleSinks)
print('The optimal flow F1 has value {:.1f}'.format(F1.value))
print('The optimal flow F2 has value {:.1f}'.format(F2.value))

---------------------
optimization problem (LP):
6 variables, 12 affine constraints

F1      : (1, 1), continuous
F2      : (1, 1), continuous
f       : dict of 4 variables, (1, 1), continuous

          maximize F1 + F2
such that
  Feasible S-(T1, T2)-flow in G has value (F1, F2).

---------------------
The optimal flow F1 has value 2.0
The optimal flow F2 has value 2.0

A similar syntax can be used for multiple sources-single sink flows.

Quadratic constraints¶

Quadratic inequalities are entered in the following way:

>>> t**2 > 2*t - alpha + x[1]*x[2]
<Quadratic Constraint: -t² + 2·t - alpha + x[1]·x[2] ≤ 0>
>>> (t & alpha) * A[1] * x + (x +2 | Z[1][:,1]) < 3*(1|Y)-alpha
<Quadratic Constraint: [t, alpha]·A[1]·x + ⟨x + [2], Z[1][:,1]⟩ - (3·⟨[1], Y⟩ - alpha) ≤ 0>

Note that PICOS does not check the convexity of convex constraints. The solver will raise an Exception if it does not support non-convex quadratics.

Second Order Cone Constraints¶

There are two types of second order cone constraints supported in PICOS.

The constraints of the type $\Vert x \Vert \leq t$ , where $t$ is a scalar affine expression and $x$ is a multidimensional affine expression (possibly a matrix, in which case the norm is Frobenius). This inequality forces the vector $[t; x]$ to belong to a Lorrentz-Cone (also called ice-cream cone).

The constraints of the type $\Vert x \Vert^2 \leq t u,\ t \geq 0$ , where $t$ and $u$ are scalar affine expressions and $x$ is a multidimensional affine expression, which constrain the vector $[t; u; x]$ inside a rotated version of the Lorretz cone.

A few examples:

>>> # A simple ice-cream cone constraint
>>> abs(x) < (2|x-1)
<5×1 SOC Constraint: ‖x‖ ≤ ⟨[2], x - [1]⟩>
>>> # SOC constraint with Frobenius norm
>>> abs(Y+Z[0].T) < t+alpha
<9×1 SOC Constraint: ‖Y + Z[0]ᵀ‖ ≤ t + alpha>
>>> # Rotated SOC constraint
>>> abs(Z[1][:,0])**2 < (2*t-alpha)*(x[2]-x[-1])
<6×1 RSOC Constraint: ‖Z[1][:,0]‖² ≤ (2·t - alpha)·(x[2] - x[-1])>
>>> # t**2 is understood as the squared norm of [t]
>>> t**2 < D['Elisa']+t
<3×1 RSOC Constraint: ‖t‖² ≤ D[Elisa] + t>
>>> # 1 is understood as the squared norm of [1]
>>> 1 < (t-1)*(x[2]+x[3])
<3×1 RSOC Constraint: ‖1^(1/2)‖² ≤ (t - 1)·(x[2] + x[3])>

Semidefinite Constraints¶

Linear matrix inequalities (LMI) can be entered thanks to an overload of the operators << and >>. For example, the LMI

$\begin{equation*} \sum_{i=0}^3 x_i b_i b_i^T \succeq b_4 b_4^T, \end{equation*}$

where $\succeq$ is used to denote the Löwner ordering, is passed to PICOS by writing:

>>> pic.sum([x[i]*b[i]*b[i].T for i in range(4)]) >> b[4]*b[4].T
<4×4 LMI Constraint: ∑(x[i]·b[i]·b[i]ᵀ : i ∈ [0…3]) ≽ b[4]·b[4]ᵀ>

Note the difference with

>>> pic.sum([x[i]*b[i]*b[i].T for i in range(4)]) > b[4]*b[4].T
<4×4 Affine Constraint: ∑(x[i]·b[i]·b[i]ᵀ : i ∈ [0…3]) ≥ b[4]·b[4]ᵀ>

which yields an elementwise inequality.

For convenience, it is possible to add a symmetric matrix variable X, by specifying the option vtype=symmetric. This has the effect to store all the affine expressions which depend on X as a function of its lower triangular elements only.

>>> sdp = pic.Problem()
>>> X = sdp.add_variable('X',(4,4),vtype='symmetric')
>>> C = sdp.add_constraint(X >> 0)
>>> print(sdp)   #doctest: +NORMALIZE_WHITESPACE
---------------------
optimization problem (SDP):
10 variables, 0 affine constraints, 10 vars in 1 SD cones
<BLANKLINE>
X   : (4, 4), symmetric
<BLANKLINE>
    find vars
such that
  X ≽ 0
---------------------

In this example, you see indeed that the problem has 10=(4*5)/2 variables, which correspond to the lower triangular elements of X.

Warning

When a constraint of the form A >> B is passed to PICOS, it is not assumed that A-B is symmetric. Instead, the symmetric matrix whose lower triangular elements are those of A-B is forced to be positive semidefnite. So, in the cases where A-B is not implicitely forced to be symmetric, you should add a constraint of the form A-B==(A-B).T in the problem.

Inequalities involving geometric means¶

It is possible to enter an inequality of the form

$t \leq \prod_{i=1}^n x_i^{1/n}$

in PICOS, where $t$ is a scalar affine expression and $x$ is an affine expression of dimension $n$ (possibly a matrix, in which case $x_i$ is counted in column major order). This inequality is internally converted to an equivalent set of second order cone inequalities, by using standard techniques (cf. e.g. [1]).

Many convex constraints can be formulated using inequalities that involve a geometric mean. For example, $t \leq x_1^{2/3}$ is equivalent to $t \leq t^{1/4} x_1^{1/4} x_1^{1/4}$ , which can be entered in PICOS thanks to the function picos.geomean() :

>>> t < pic.geomean(t //x[1] //x[1] //1)
<Geometric Mean Constraint: t ≤ geomean([t; x[1]; x[1]; 1])>

Note that the latter example can also be passed to picos in a more simple way, thanks to an overloading of the ** exponentiation operator:

>>> t < x[1]**(2./3)
<Power Constraint: x[1]^(2/3) ≥ t>

Inequalities involving geometric means are stored in a temporary object of the class GeoMeanConstraint, which can be passed to a problem with add_constraint():

>>> geom_ineq = prob.add_constraint(t<pic.geomean(Y[:6]))

The constraint has an attribute constraints which contains auxiliary SOC inequalities that are used internally to represent the geometric mean:

>>> pprint(geom_ineq.constraints)
(<3×1 RSOC Constraint: ‖u[1:0-1]‖² ≤ Y[:6][0]·Y[:6][1]>,
 <3×1 RSOC Constraint: ‖u[1:2-3]‖² ≤ Y[:6][2]·Y[:6][3]>,
 <3×1 RSOC Constraint: ‖u[1:4-5]‖² ≤ Y[:6][4]·Y[:6][5]>,
 <3×1 RSOC Constraint: ‖u[2:0-3]‖² ≤ u[1:0-1]·u[1:2-3]>,
 <3×1 RSOC Constraint: ‖u[2:4-x]‖² ≤ u[1:4-5]·t>,
 <3×1 RSOC Constraint: ‖t‖² ≤ u[2:0-3]·u[2:4-x]>)

Likewise, auxiliary variables can be accessed via the variables attribute:

>>> pprint(geom_ineq.variables)
(<1×1 Continuous Variable: u[1:0-1]>,
 <1×1 Continuous Variable: u[1:2-3]>,
 <1×1 Continuous Variable: u[1:4-5]>,
 <1×1 Continuous Variable: u[2:0-3]>,
 <1×1 Continuous Variable: u[2:4-x]>)

Inequalities involving real powers or trace of matrix powers¶

As mentionned above, the ** exponentiation operator has been overloaded to support real exponents. A rational approximation of the exponent is used, and the inequality are internally reformulated as a set of equivalent SOC inequalities. Note that only inequalities defining a convex regions can be passed:

>>> t**0.6666 > x[0]
<Power Constraint: t^(2/3) ≥ x[0]>
>>> t**-0.5 < x[0]
<Power Constraint: t^(-1/2) ≤ x[0]>
>>> try:
...     t**-0.5 > x[0]
... except Exception as ex:
...     print('Exception: '+str(ex)) #doctest: +NORMALIZE_WHITESPACE
Exception: May only lower bound a concave power (p >= 0 and p <= 1).
>>> t**2 < x[1]+x[2]
<3×1 RSOC Constraint: ‖t‖² ≤ x[1] + x[2]>

More generally, inequalities involving trace of matrix powers can be passed to PICOS, by using the picos.tracepow() function. The following example creates the constraint

$\operatorname{trace}\ \big(x_0 A_0 A_0^T + x_2 A_2 A_2^T\big)^{2.5} \leq 3.$

>>> pic.tracepow(x[0] * A[0]*A[0].T + x[2] * A[2]*A[2].T, 2.5) <= 3
<Trace of Power Constraint: trace((x[0]·A[0]·A[0]ᵀ + x[2]·A[2]·A[2]ᵀ)^(5/2)) ≤ 3>

Warning

when a power expression $x^p$ (resp. the trace of matrix power $\operatorname{trace}\ X^p$ ) is used, the base $x$ is forced to be nonnegative (resp. the base $X$ is forced to be positive semidefinite) by picos.

When the exponent is $0<p<1$ , it is also possible to represent constraints of the form $\operatorname{trace}(M X^p) \geq t$ with SDPs, where $M\succeq 0$ , see [2].

>>> pic.tracepow(X, 0.6666, coef = A[0].T*A[0]) >= t
<Trace of Power Constraint: trace(A[0]ᵀ·A[0]·X^(2/3)) ≥ t>

As for geometric means, inequalities involving real powers yield their internal representation via the constraints and variables attributes.

Inequalities involving generalized p-norm¶

Inequalities of the form $\Vert x \Vert_p \leq t$ can be entered by using the function picos.norm(). This function is also defined for $p < 1$ by the usual formula $\Vert x \Vert_p := \Big(\sum_i |x_i|^p \Big)^{1/p}$ . The norm function is convex over $\mathbb{R}^n$ for all $p\geq 1$ , and concave over the set of vectors with nonnegative coordinates for $p \leq 1$ .

>>> pic.norm(x,3) < t
<p-Norm Constraint: ‖x‖_3 ≤ t>
>>> pic.norm(x,'inf') < 2
<p-Norm Constraint: ‖x‖_inf ≤ 2>
>>> pic.norm(x,0.5) > x[0]-x[1]
<Generalized p-Norm Constraint: ‖x‖_(1/2) ≥ x[0] - x[1]>

Warning

Note that when a constraint of the form norm(x,p) >= t is entered (with $p \leq 1$ ), PICOS forces the vector x to be nonnegative (componentwise).

Inequalities involving the generalized $L_{p,q}$ norm of a matrix can also be handled with picos, cf. the documentation of picos.norm() .

As for geometric means, inequalities involving p-norms yield their internal representation via the constraints and variables attributes.

Inequalities involving the nth root of a determinant¶

The function picos.detrootn() can be used to enter the $n$ -th root of the determinant of a $(n \times n)$ -symmetric positive semidefinite matrix:

>>> M = sdp.add_variable('M',(5,5),'symmetric')
>>> t < pic.detrootn(M)
<n-th Root of a Determinant Constraint: det(M)^(1/5) ≥ t>

Warning

Note that when a constraint of the form t < pic.detrootn(M) is entered (with $p \leq 1$ ), PICOS forces the matrix M to be positive semidefinite.

As for geometric means, inequalities involving the nth root of a determinant yield their internal representation via the constraints and variables attributes.

Set membership¶

Since Picos 1.0.2, there is a Set class that can be used to pass constraints as membership of an affine expression to a set.

Following sets are currently supported:

$L_p-$ balls representing the set $\{x: \Vert x \Vert_p \leq r \}$ can be constructed with the function pic.ball()

The standard simplex (scaled by a factor $\gamma$ ) $\{x \geq 0: \sum_i x_i \leq r \}$ can be constructed with the function pic.simplex()

Truncated simplexes $\{0 \leq x \leq 1: \sum_i x_i \leq r \}$ and symmetrized Truncated simplexes $\{x: \Vert x \Vert_\infty \leq 1, \Vert x \Vert_1\leq r \}$ can be constructed with the function pic.truncated_simplex()

Membership of an affine expression to a set can be expressed with the overloaded operator <<. This returns a temporary object that can be passed to a picos problem with the function add_constraint() .

>>> x << pic.simplex(1)
<Standard Simplex Constraint: x ∈ {x ≥ 0 : ∑(x) ≤ 1}>
>>> x << pic.truncated_simplex(2)
<Truncated Simplex Constraint: x ∈ {0 ≤ x ≤ 1 : ∑(x) ≤ 2}>
>>> x << pic.truncated_simplex(2,sym=True)
<Symmetrized Truncated Simplex Constraint: x ∈ {-1 ≤ x ≤ 1 : ∑(|x|) ≤ 2}>
>>> x << pic.ball(3)
<5×1 SOC Constraint: ‖x‖ ≤ 3>
>>> pic.ball(2,'inf') >> x
<p-Norm Constraint: ‖x‖_inf ≤ 2>
>>> x << pic.ball(4,1.5)
<p-Norm Constraint: ‖x‖_(3/2) ≤ 4>

Write a Problem to a file¶

It is possible to write a problem to a file, thanks to the function write_to_file(). Several file formats and file writers are available, have a look at the doc of write_to_file() for more explanations.

Below is a hello world example, which writes a simple LP to a .lp file:

import picos as pic
prob = pic.Problem(verbose = False)
y = prob.add_variable('y',1)
x = prob.add_variable('x',1)
prob.add_constraint(x>1.5)
prob.add_constraint(y-x>0.7)
prob.set_objective('min',y)
#let first picos display the problem
print(prob)
print()
#now write the problem to a .lp file...
prob.write_to_file('helloworld.lp')
print()
#and display the content of the freshly created file:
print(open('helloworld.lp').read())

Generated output:

---------------------
optimization problem (LP):
2 variables, 2 affine constraints

x   : (1, 1), continuous
y   : (1, 1), continuous

        minimize y
such that
x ≥ 1.5
y - x ≥ 0.7
---------------------

writing problem in helloworld.lp...
done.

\* file helloworld.lp generated by picos*\
Minimize
obj : 1 y
Subject To
in0 : -1 x <= -1.5
in1 : -1 y+ 1 x <= -0.7
Bounds
y free
x free
Generals
Binaries
End

Solve a Problem¶

To solve a problem, you have to use the method solve() of the class Problem. Alternatively, the functions maximize(obj) and minimize(obj) can be used to specify the objective function and call the solver in a single statement. These method accept several options. In particular the solver can be specified by passing an option of the form solver='solver_name'. For a list of available parameters with their default values, see the doc of the function set_all_options_to_default().

Once a problem has been solved, the optimal values of the variables are accessible with the value property. Depending on the solver, you can also obtain the slack and the optimal dual variables of the constraints thanks to the properties dual and slack of the class Constraint. See the doc of dual for more explanations on the dual variables for second order cone programs (SOCP) and semidefinite programs (SDP).

The class Problem also has two interesting properties: type, which indicates the class of the optimization problem (‘LP’, ‘SOCP’, ‘MIP’, ‘SDP’,…), and status, which indicates if the problem has been solved (the default is 'unsolved'; after a call to solve() this property can take the value of any code returned by a solver, such as 'optimal', 'unbounded', 'near-optimal', 'primal infeasible', 'unknown', …).

Below is a simple example, to solve the linear programm:

$\begin{center} $\begin{array}{ccc} \underset{x \in \mathbb{R}^2}{\mbox{minimize}} & 0.5 x_1 + x_2 &\\ \mbox{subject to} & x_1 &\geq x_2\\ & \left[ \begin{array}{cc} 1 & 0\\ 1 & 1 \end{array} \right] x &\leq \left[ \begin{array}{c} 3 \\4 \end{array} \right]. \end{array}$ \end{center}$

More examples can be found here.

P = pic.Problem()
A = pic.new_param('A', cvx.matrix([[1,1],[0,1]]) )
x = P.add_variable('x',2)
P.add_constraint(x[0]>x[1])
P.add_constraint(A*x<[3,4])
objective = 0.5 * x[0] + x[1]
P.set_objective('max', objective) #or directly P.maximize(objective)

#display the problem and solve it
print(P)
print('type:   '+P.type)
print('status: '+P.status)
P.solve(verbose=0,solver='cvxopt')
print('status: '+P.status)
print()

#--------------------#
#  objective value   #
#--------------------#

print('The optimal value of this problem is:')
print(P.obj_value()) # "print(objective)" would also work, because objective is valued
print()

#--------------------#
#  optimal variable  #
#--------------------#
x_opt = x.value
print('The solution of the problem is:')
print(x_opt) # "print(x)" would also work, since x is now valued

#--------------------#
#  slacks and duals  #
#--------------------#
c0=P.get_constraint(0)
print('The dual of the constraint {} is:'.format(c0))
print(c0.dual)
print()
print('And its slack is:')
print(c0.slack)
print()

c1=P.get_constraint(1)
print('The dual of the constraint {} is:'.format(c1))
print(c1.dual)
print('And its slack is:')
print(c1.slack)

---------------------
optimization problem (LP):
2 variables, 3 affine constraints

x   : (2, 1), continuous

    maximize 0.5·x[0] + x[1]
such that
  x[0] ≥ x[1]
  A·x ≤ [2×1]
---------------------
type:   LP
status: unsolved
status: optimal

The optimal value of this problem is:
3.0...

The solution of the problem is:
[ 2.00e+00]
[ 2.00e+00]

The dual of the constraint x[0] ≥ x[1] is:
0.25...

And its slack is:
1.8...

The dual of the constraint A·x ≤ [2×1] is:
[ 4.56e-10]
[ 7.50e-01]

And its slack is:
[ 1.00e+00]
[-8.71e-10]

A note on dual variables¶

For second order cone constraints of the form $\Vert \mathbf{x} \Vert \leq t$ , where $\mathbf{x}$ is a vector of dimension $n$ , the dual variable is a vector of dimension $n+1$ of the form $[\lambda; \mathbf{z}]$ , where the $n$ -dimensional vector $\mathbf{z}$ satisfies $\Vert \mathbf{z} \Vert \leq \lambda$ .

For rotated second order cone constraints of the form $\Vert \mathbf{x} \Vert^2 \leq t u$ where $\mathbf{x}$ is a vector of dimension $n$ , the dual variable is a vector of dimension $n+2$ of the form $[\alpha, \beta, \mathbf{z}]$ with $\alpha, \beta \geq 0$ , where the $n$ -dimensional vector $\mathbf{z}$ satisfies $\Vert \mathbf{z} \Vert^2 \leq 4 \alpha \beta$ .

In general, a linear problem with second order cone constraints (both standard and rotated) and semidefinite constraints can be written under the form:

$\begin{center} $\begin{array}{cclc} \underset{\mathbf{x} \in \mathbb{R}^n}{\mbox{minimize}} & \mathbf{c}^T \mathbf{x} & &\\ \mbox{subject to} & A^e \mathbf{x} + \mathbf{b^e} & = 0 &\\ & A^l \mathbf{x} + \mathbf{b^l} & \leq 0 &\\ & \Vert A^s_i \mathbf{x} + \mathbf{b^s_i} \Vert &\leq \mathbf{f^s_i}^T \mathbf{x} +d^s_i, & \forall i \in I\\ & \Vert A^r_j \mathbf{x} + \mathbf{b^r_j} \Vert^2 &\leq (\mathbf{f^{r_1}_j}^T \mathbf{x} +d^{r_1}_j) (\mathbf{f^{r_2}_j}^T \mathbf{x} +d^{r_2}_j), & \forall j \in J\\ & 0 & \leq \mathbf{f^{r_1}_j}^T \mathbf{x} +d^{r_1}_j, & \forall j \in J\\ & \sum_{i=1}^n x_i M_i & \succeq M_0 \end{array}$ \end{center}$

where

$\mathbf{c}, \big\{\mathbf{f^s_i}\big\}_{i\in I}, \big\{\mathbf{f^{r_1}_j}\big\}_{j \in J}, \big\{\mathbf{f^{r_2}_j}\big\}_{j \in J}$ are vectors of dimension $n$ ;

$\big\{d^s_i\big\}_{i \in I}, \big\{d^{r_1}_j\big\}_{j \in J}, \big\{d^{r_2}_j\big\}_{j \in J}$ are scalars;

$\big\{\mathbf{b^s_i}\big\}_{i\in I}$ are vectors of dimension $n^s_i$ and $\big\{A^s_i\big\}_{i\in I}$ are matrices of size $n^s_i \times n$ ;

$\big\{\mathbf{b^r_j}\big\}_{j\in J}$ are vectors of dimension $n^r_j$ and $\big\{A^r_j\big\}_{j\in J}$ are matrices of size $n^r_j \times n$ ;

$\mathbf{b^e}$ is a vector of dimension $n^e$ and $A^e$ is a matrix of size $n^e \times n$ ;

$\mathbf{b^l}$ is a vector of dimension $n^l$ and $A^l$ is a matrix of size $n^l \times n$ ;

$\big\{M_k\big\}_{k=0,\ldots,n}$ are $m \times m$ symmetric matrices ( $M_k \in \mathbb{S}_m$ ).

Its dual problem can be written as:

$\begin{center} $\begin{array}{cll} \mbox{maximize} & \mathbf{b^e}^T \mathbf{\mu^e} + \mathbf{b^l}^T \mathbf{\mu^l} & + \sum_{i\in I} \big( \mathbf{b^s_i}^T \mathbf{z^s_i} - d^s_i \lambda_i \big) + \sum_{j\in J} \big( \mathbf{b^r_j}^T \mathbf{z^r_j} - d^{r_1}_j \alpha_j - d^{r_2}_j \beta_j \big) + \langle M_0, X \rangle\\ \mbox{subject to} & c + {A^e}^T \mathbf{\mu^e} + {A^l}^T \mathbf{\mu^l} & + \sum_{i\in I} \big( {A^s_i}^T \mathbf{z^s_i} -\lambda_i \mathbf{f^s_i} \big) + \sum_{j\in J} \big( {A^r_j}^T \mathbf{z^r_j} -\alpha_j \mathbf{f^{r_1}_j} - \beta_j \mathbf{f^{r_2}_j} \big) = \mathcal{M} \bullet X \\ & \mu_l \geq 0 &\\ & \Vert \mathbf{z^s_i} \Vert \leq \lambda_i, &\forall i \in I\\ & \Vert \mathbf{z^r_j} \Vert^2 \leq 4 \alpha_j \beta_j, &\forall j \in J\\ & \ \ 0 \ \ \ \leq \alpha_j, &\forall j \in J\\ & X \succeq 0 \end{array}$ \end{center}$

where $\mathcal{M} \bullet X$ stands for the vector of dimension $n$ with $\langle M_i, X \rangle$ on the $i$ th coordinate, and the dual variables are

$\mu^e \in \mathbb{R}^{n_e}$

$\mu^l \in \mathbb{R}^{n_l}$

$z^s_i \in \mathbb{R}^{n^s_i},\ \forall i \in I$

$\lambda_i \in \mathbb{R},\ \forall i \in I$

$z^r_j \in \mathbb{R}^{n^r_j},\ \forall j \in J$

$(\alpha_j,\beta_j) \in \mathbb{R}^2,\ \forall j \in J$

$X \in \mathbb{S}_m$

When quering the dual of a constraint of the above primal problem, picos will return

$\mu^e$ for the constraint $A^e \mathbf{x} + \mathbf{b^e} = 0$ ;

$\mu^l$ for the constraint $A^l \mathbf{x} + \mathbf{b^l} \geq 0$ ;

The $(n^s_i+1)$ -dimensional vector $\mu^s_i = [\lambda_i;\mathbf{z^s_i}]$ for the constraint

$\Vert A^s_i \mathbf{x} + \mathbf{b^s_i} \Vert \leq \mathbf{f^s_i}^T \mathbf{x} +d^s_i$ ;

The $(n^r_j+2)$ -dimensional vector $\mu^r_j = [\alpha_j; \beta_j; \mathbf{z^r_j}]$ for the constraint

$\Vert A^r_j \mathbf{x} + \mathbf{b^r_j} \Vert^2 \leq (\mathbf{f^{r_1}_j}^T \mathbf{x} +d^{r_1}_j) (\mathbf{f^{r_2}_j}^T \mathbf{x} +d^{r_2}_j)$ ;

The symmetric positive definite matrix $X$ for the constraint

$\sum_{i=1}^n x_i M_i \succeq M_0$ .

References¶

“Applications of second-order cone programming”, M.S. Lobo, L. Vandenberghe, S. Boyd and H. Lebret, Linear Algebra and its Applications, 284, p. 193-228, 1998.

“On the semidefinite representations of real functions applied to symmetric matrices”, G. Sagnol, Linear Algebra and its Applications, 439(10), p. 2829-2843, 2013.