Firstly, I will present the problem.

Maximize vector function with constraints:

We have vector function. $\phi$ and $p$ are n-dimensional vectors. With given n (or if it is not possible - with n=2), find real values, that will maximize globally this function.

$$ f(\phi, p, \theta) = 1-\sum_{i=1}^{n} w_i (-\lambda_i \log_2\lambda_i - (1-\lambda_i) log_2(1-\lambda_i)) $$

$$ w_i = (1+\cos\phi_i \cos\theta)p_i $$

$$ \lambda_i = \frac{ 1 + cos(\phi_i - \theta) }{2(1+cos\phi_i cos\theta)} $$


$$ \sum_{i=1}^{n} p_i =1 $$

$$ \forall i\ldotp p_i >0 $$

$$ \sum_{i=1}^{n} w_i \lambda_i =\frac{1}{2} $$

$$ \forall i\ldotp \lambda_i >0 $$

Of course we should also take care about constraint

$$ \forall i\ldotp cos\phi_i cos\theta \neq -1 $$

to prevent the situation, when denominator is 0.

My attempt

I tried to solve this problem with scipy.optimize package. I had to pass the initial vector (initial guess). And value, which was given by algorithm was dependent on this initial guess. So I tried to take random initial guesses and tried to find maximal values among more and more initial guesses

import numpy as np
from scipy.optimize import minimize

def J(data):
    l = len(data)
    if (l%2==0 | l<2):
        raise Exception('Wrong length or arguments!')

    phi = data[0:l//2]
    p = data[l//2:l-1]
    theta = data[l-1]

    phi = np.array(phi)    
    p = np.array(p)

    w = (1 + np.cos(phi)*np.cos(theta))*p
    lam = (1+np.cos(phi-theta))/(2*(1+np.cos(phi)*np.cos(theta)))

    return -(1 - sum(w * ((-lam*np.log2(lam)) - (1-lam)*np.log2(1-lam))))

def probability_constraint(x):
    return sum(x[len(x)//2:len(x)-1])-1

def lambda_w_constraint(data):
    l = len(data)
    phi = data[0:l//2]
    p = data[l//2:-1]
    theta = data[l-1]

    phi = np.array(phi)    
    p = np.array(p)

    w = (1 + np.cos(phi)*np.cos(theta))*p
    lam = (1+np.cos(phi-theta))/(2*(1+np.cos(phi)*np.cos(theta)))    
    return sum(w*lam)-1/2

def lambda_greater_than_zero(data):
    l = len(data)
    phi = data[0:l//2]
    p = data[l//2:-1]
    theta = data[l-1]

    lam = (1+np.cos(phi-theta))/(2*(1+np.cos(phi)*np.cos(theta)))    
    return lam

def probability_positive(x):
    return x

cons = ({'type': 'ineq',
'fun' : lambda x: lambda_greater_than_zero(x)},
        {'type': 'eq',
'fun' : lambda x: probability_constraint(x)},
        {'type': 'eq',
'fun' : lambda x: lambda_w_constraint(x)},
'fun' : lambda x: probability_positive(x)}


def maximize(numiter):
    maxJ = -1000
    maxx = []

    for i in range(0, numiter):
        x0 = np.random.rand(5)*10
        res = minimize(J, x0, constraints=cons)
        r = J(res.x)
        if (r>maxJ):
            maxJ = r
            maxx = res.x

    return maxJ, maxx



I know, that this method is very naive and incorrect.

So my question is: do you know how to solve this problem? It can be in python3, but other languages (Matlab?) are also good for me.


Tried it in MATLAB with the modelling toolbox YALMIP (disclaimer, developed by me)

theta = sdpvar(1);
n = 5;
p = sdpvar(n,1);
f = sdpvar(n,1);
w = sdpvar(n,1);
l = sdpvar(n,1);

Model = [w == (1 + cos(f)*cos(theta)).*p,
         l == (1 + cos(f-theta))./(2 + 2*cos(f)*cos(theta)),        
    sum(p) == 1, p>=0, sum(w.*p) == .5, l >= 0]
Model = [Model, -pi <= theta <= pi, -pi <= f <= pi];
objective = 0;
for i = 1:n
   objective = objective+w(i)*(entropy(l(i)) + entropy(1-l(i)));     
optimize(Model, objective)

Both solvers fmincon and ipopt returned reasonable solutions. I use the entropy operator $entropy(x) = -x\log(x)$ as YALMIP can exploit some properties better when computing the derivatives and avoids some issues around $0$.

You can also test the global solver available in YALMIP

optimize(Model, objective, sdpsettings('solver','bmibnb'))

At least it appears to find the global solution for $n=2$.


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.