0

I'm trying to optimize function objective.

full_constr_data consists of 6 type of goals, each goal is divided by years, each year is represented by project-based data. So, I'm weighting full_constr_data project-wise by argument x of function objective, e.g. full_constr_data[0][2][3] * x[3] means that "goal #0, year #2, project #3 is weighted by x[3] The results are stored in variable full_constr_data_weighted.

The next step is project-wise sum up of full_constr_data_weighted. Example, sum all projects for goal #0, year #2:

full_constr_data_weighted[0][2][0] + full_constr_data_weighted[0][2][1] + ... + full_constr_data_weighted[0][2][n]

where 'n' - total number of projects. The data is stored in variable full_sum.

After that I'm calculating probabilities. I take quantiles from variable constr_mod and based on the value calculate the probabilities of exceeding the quantile. constr_mod and full_sum have exactly the same structure, however for each goal # and following year # constr_mode contains a single value, while full_sum has vector of values (distribution). Calculated probabilities are stored to variable my_prob.

Finally, I'm summing up all the probabilities in my_prob. This sum I have to optimize: make it as large as possible (notice a minus sign in return statement).

Optimization problem has a single inequality constraint: sum of vector obj weighted by x should be larger that 1000. Interpret obj as NPV for each of the project.

I use diffev2 from package mystic.

Variables full_constr_data, constr_mod, pen_mult, obj are stored in my_data.spydata file: download via Google Drive

Unfortunately, optimization didn't converge:

Generation 11980 has ChiSquare: inf
Generation 11990 has ChiSquare: inf
Generation 12000 has ChiSquare: inf
STOP("EvaluationLimits with {'evaluations': 1200000, 'generations': 12000}")

Any suggestions how to solve this non-convex problem?

import numpy as np
from mystic.solvers import diffev2
from mystic.monitors import VerboseMonitor
import mystic.symbolic as ms

def objective(x):
    # 'full_constr_data' weighted by argument 'x'
    full_constr_data_weighted = []
    for i in range(len(full_constr_data)):
        temp = []
        for k in range(len(full_constr_data[i])):
            temp.append( [ full_constr_data[i][k][m] * x[m] \
                          for m in range(len(full_constr_data[i][k])) ] )
        full_constr_data_weighted.append(temp)

    # Project-wise sum of weighted data
    full_sum = []
    for i in range(len(full_constr_data_weighted)):
        temp = []
        for j in range(len(full_constr_data_weighted[i])):
            temp2 = np.array( [ 0 for m in range(len(full_constr_data_weighted[i][j][0])) ] )
            for k in range(len(full_constr_data_weighted[i][j])):
                temp2 = temp2 + full_constr_data_weighted[i][j][k]
            temp.append(temp2)
        full_sum.append(temp)

    # Probability calculation
    my_prob = []
    for i in range(len(full_sum)):
        temp = []
        for j in range(len(full_sum[i])):
            temp.append(sum(full_sum[i][j] > constr_mod[i][j]) / len(full_sum[i][j]))
        my_prob.append(np.array(temp))

    # Probability data weighted by 'pen_mult'
    my_prob_uweighted = list(np.array(my_prob) * np.array(pen_mult))

    # Sum of all weighted probability data (function to maximize)
    sum_prob = sum([sum(my_prob_uweighted[i]) for i in range(len(my_prob_uweighted))])

    return -sum_prob



# Inequality constraint
equation = 'x0*{0} + x1*{1} + x2*{2} + x3*{3} + x4*{4} + x5*{5} + x6*{6} + x7*{7} + x8*{8} + x9*{9} + x10*{10} + x11*{11} + x12*{12} + x13*{13} + x14*{14} + x15*{15} + x16*{16} + x17*{17} + x18*{18} + x19*{19} + x20*{20} + x21*{21} + x22*{22} + x23*{23} + x24*{24} + x25*{25} + x26*{26} + x27*{27} + x28*{28} + x29*{29} >= {30}'\
.format(obj[0], obj[1], obj[2], obj[3], obj[4], obj[5], obj[6], obj[7], obj[8], obj[9],
        obj[10], obj[11], obj[12], obj[13], obj[14], obj[15], obj[16], obj[17], obj[18], obj[19],
        obj[20], obj[21], obj[22], obj[23], obj[24], obj[25], obj[26], obj[27], obj[28], obj[29],
        1000)


cf = ms.generate_constraint(ms.generate_solvers(ms.simplify(equation)))

bounds = [(0,1)]*30

mon = VerboseMonitor(10)

result = diffev2(objective, x0=bounds, bounds=bounds, constraints=cf, \
                 npop=40, gtol=200, disp=False, full_output=True, itermon=mon)
kimsergeo
  • 1
  • 1
  • Hi, I'm the mystic author. I didn't run your code, since to do so, I'd have to download stuff from the google drive. However, if you are seeing `inf` everywhere, then that typically means that you have constraints that have no solution. Try passing a test solution vector directly to the constraints function. Alternately, you can build a penalty instead of a constraint. Then see what happens. I'll see if I can mock your problem up with dummy data in the meantime. – Mike McKerns May 04 '19 at 11:44
  • Again, it's not that the solver didn't converge, it actually didn't see any solutions to the constraints. With a dummy objective and data, your constraints work as expected and find a solution. I tried checking `cf(np.random.rand(30))`, and the constraints hold. If you have a particularly difficult nonlinear surface in the objective, I'd also try one of the `ensemble` solvers, with a `Powell` solver being run in parallel. Or try with a `penalty`, as I suggested earlier. The point being, `inf` means no valid points found (yet), due to the bounds, constraints, and objective. – Mike McKerns May 04 '19 at 13:21
  • If you post minimal working code, that shows the same behavior, you'll likely get more help. – Mike McKerns May 04 '19 at 13:22
  • I tried to use ```penalty``` instead of inequality constraints and faced strange behavior. Negative ```ChiSquare``` values were displayed, and they continue to reduce (e.g. ```-12```, ```-34```, etc.). As far as I know Chi-Square distribution has support of positive real numbers. Is it solver's mistake, or maybe I interpreted ```ChiSquare``` wrongly? – kimsergeo May 09 '19 at 14:21
  • If you don't modify what the monitor saves, then the monitor is saving `x` (solution vector) versus `y` (output of the cost function). I'm assuming the negative values you are experiencing are simply the penalized negative cost. – Mike McKerns May 09 '19 at 17:09

0 Answers0