CODEWITHSUNDEEP
optimizationvectorgeometryconstraintsmathematical-optimization
gokva
gokva

Reputation: 55

Formulating a constraint on the SUM of the magnitudes of a set of vectors in a constrained optimization

I'm solving a geometric constrained optimization problem. The variables in the optimization are the x-y components for a set of vectors. The objective function is quadratic in these variables.

However, I need to constrain the SUM of the magnitudes of a subset of the vectors.

Specifically, suppose this subset consists of {v1,v2,...,vn}

I need the solution to satisfy

||v1|| + ||v2|| + .... + ||vn|| < L

If it was just a single vector I could square both sides to get a quadratic constraint and frame the problem as a QCQP

v1.x * v1.x + v1.y * v1.y < L*L

However, I have multiple vectors. So is there any way to express the constraint in such a way that I could apply a technique more specific than general non-linear constrained optimization? Or given that my objective function can be minimized analytically, would it make sense to solve the problem by

  1. Ignoring the constraint and obtaining an optimum value x* for the objective function
  2. Projecting x* onto the constraint manifold numerically to get a final solution that satisfies the constraints?

Upvotes: 0

Views: 784

Answers (1)

sascha
sascha

Reputation: 33522

Not sure what else is within your optimization-problem, but the task of constraining your norms alone is nonlinear but convex which allows efficient solving.

Using external libraries you can prototype it like this. Here cvxpy (python) is used.

There are many similar libraries following the same ideas like: cvxopt (python), picos (python), yalmip (matlab), convex.jl (julia). The official solver-APIs are usually way more low-level and there is more work to do. Inbetween there is also JuMP (julia). 

from cvxpy import *

L = 10.0
V = Variable(3,5)  # 3 vectors

constraints = []
constraints.append(sum_entries(norm(V, axis=1)) <= L)

objective = Maximize(sum_entries(V))
prob = Problem(objective, constraints)
prob.solve()
print("status:", prob.status)
print("optimal value", prob.value)
print("optimal var", V.value)
print('constr: ', sum_entries(norm(V, axis=1).value))

Output:

status: optimal
optimal value 22.36067971461066
optimal var [[ 1.49071198  1.49071198  1.49071198  1.49071198  1.49071198]
 [ 1.49071198  1.49071198  1.49071198  1.49071198  1.49071198]
 [ 1.49071198  1.49071198  1.49071198  1.49071198  1.49071198]]
constr:  sum_entries([[ 3.33333332]
 [ 3.33333332]
 [ 3.33333332]])

The above is automatically converted so SOCP-form and can be solved by commercial-solvers or open-source solvers like ECOS and SCS.

This conversion also proves to us, that this problem is convex (by construction)! The approach is called Disciplined Convex Programming .

Depending on your library/software chosings, you have to do this conversion manually. It should be not that hard when you introduce some helper-variables to collect the norms of your vectors. Within Gurobi, you just would need to use the basic SOCP-constraint docs.

Remark: ||v1|| + ||v2|| + .... + ||vn|| < L is scary, as numerical-optimization usually only cares about <=. Everything else needs trickery (epsilon-values...)

Edit:

Here is a pure Gurobi-approach, which can give you some idea on how to achieve this with more low-level libraries supporting similar functions as Gurobi's API (i'm thinking about Mosek and CPLEX without knowing there APIs much; i think Mosek is quite different).

from gurobipy import *
import numpy as np

L = 10.0

# Model
m = Model("test")

# Vars
v = np.empty((3,5), dtype=object)
for i in range(3):
    for j in range(5):
        v[i, j] = m.addVar()  # by default >= 0; it's just an example

norm_v = np.empty(3, dtype=object)
for i in range(3):
    norm_v[i] = m.addVar()    # aux-vars to collect norms

m.update()  # make vars usable for posting constraints

# Constraints
for i in range(3):
    m.addQConstr(np.dot(v[i, :], v[i, :]),
                 GRB.LESS_EQUAL, norm_v[i] * norm_v[i])  # this is the SOCP-constraint for our norm

m.addConstr(np.sum(norm_v) <= L)  # gurobi-devs would propose using quicksum

# Objective
m.setObjective(np.sum(v), GRB.MAXIMIZE)

# Solve
m.optimize()

def get_val(x):
    return x.X
get_val_func = np.vectorize(get_val)
print('optimal var: ', get_val_func(v))

Output:

Optimal objective 2.23606793e+01

optimal var:  [[ 1.49071195  1.49071195  1.49071195  1.49071195  1.49071195]
 [ 1.49071195  1.49071195  1.49071195  1.49071195  1.49071195]
 [ 1.49071195  1.49071195  1.49071195  1.49071195  1.49071195]]

Upvotes: 1

Related Questions