# lagrangian_relaxation.py¶

This example is inspired by an entry on the “adventures in optimization” blog. The sources of the article can be found here. This example solves the generalized assignment problem, with or without Lagrangian relaxation.

  1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 # -------------------------------------------------------------------------- # Source file provided under Apache License, Version 2.0, January 2004, # http://www.apache.org/licenses/ # (c) Copyright IBM Corp. 2015, 2018 # -------------------------------------------------------------------------- import json from docplex.util.environment import get_environment from docplex.mp.model import Model # ---------------------------------------------------------------------------- # Initialize the problem data # ---------------------------------------------------------------------------- B = [15, 15, 15] C = [ [ 6, 10, 1], [12, 12, 5], [15, 4, 3], [10, 3, 9], [8, 9, 5] ] A = [ [ 5, 7, 2], [14, 8, 7], [10, 6, 12], [ 8, 4, 15], [ 6, 12, 5] ] # ---------------------------------------------------------------------------- # Build the model # ---------------------------------------------------------------------------- def run_GAP_model(As, Bs, Cs, **kwargs): with Model('GAP per Wolsey -without- Lagrangian Relaxation', **kwargs) as mdl: print("#As={}, #Bs={}, #Cs={}".format(len(As), len(Bs), len(Cs))) number_of_cs = len(C) # variables x_vars = [mdl.binary_var_list(c, name=None) for c in Cs] # constraints mdl.add_constraints(mdl.sum(xv) <= 1 for xv in x_vars) mdl.add_constraints(mdl.sum(x_vars[ii][j] * As[ii][j] for ii in range(number_of_cs)) <= bs for j, bs in enumerate(Bs)) # objective total_profit = mdl.sum(mdl.scal_prod(x_i, c_i) for c_i, x_i in zip(Cs, x_vars)) mdl.maximize(total_profit) # mdl.print_information() s = mdl.solve() assert s is not None obj = s.objective_value print("* GAP with no relaxation run OK, best objective is: {:g}".format(obj)) return obj def run_GAP_model_with_Lagrangian_relaxation(As, Bs, Cs, max_iters=101, **kwargs): with Model('GAP per Wolsey -with- Lagrangian Relaxation', **kwargs) as mdl: print("#As={}, #Bs={}, #Cs={}".format(len(As), len(Bs), len(Cs))) number_of_cs = len(Cs) c_range = range(number_of_cs) # variables x_vars = [mdl.binary_var_list(c, name=None) for c in Cs] p_vars = mdl.continuous_var_list(Cs, name='p') # new for relaxation mdl.add_constraints(mdl.sum(xv) == 1 - pv for xv, pv in zip(x_vars, p_vars)) mdl.add_constraints(mdl.sum(x_vars[ii][j] * As[ii][j] for ii in c_range) <= bs for j, bs in enumerate(Bs)) # lagrangian relaxation loop eps = 1e-6 loop_count = 0 best = 0 initial_multiplier = 1 multipliers = [initial_multiplier] * len(Cs) total_profit = mdl.sum(mdl.scal_prod(x_i, c_i) for c_i, x_i in zip(Cs, x_vars)) mdl.add_kpi(total_profit, "Total profit") while loop_count <= max_iters: loop_count += 1 # rebuilt at each loop iteration total_penalty = mdl.scal_prod(p_vars, multipliers) mdl.maximize(total_profit + total_penalty) s = mdl.solve() if not s: print("*** solve fails, stopping at iteration: %d" % loop_count) break best = s.objective_value penalties = [pv.solution_value for pv in p_vars] print('%d> new lagrangian iteration:\n\t obj=%g, m=%s, p=%s' % (loop_count, best, str(multipliers), str(penalties))) do_stop = True justifier = 0 for k in c_range: penalized_violation = penalties[k] * multipliers[k] if penalized_violation >= eps: do_stop = False justifier = penalized_violation break if do_stop: print("* Lagrangian relaxation succeeds, best={:g}, penalty={:g}, #iterations={}" .format(best, total_penalty.solution_value, loop_count)) break else: # update multipliers and start loop again. scale_factor = 1.0 / float(loop_count) multipliers = [max(multipliers[i] - scale_factor * penalties[i], 0.) for i in c_range] print('{0}> -- loop continues, m={1!s}, justifier={2:g}'.format(loop_count, multipliers, justifier)) return best def run_default_GAP_model_with_lagrangian_relaxation(**kwargs): return run_GAP_model_with_Lagrangian_relaxation(As=A, Bs=B, Cs=C, **kwargs) # ---------------------------------------------------------------------------- # Solve the model and display the result # ---------------------------------------------------------------------------- if __name__ == '__main__': # Run the model. If a key has been specified above, the model will run on # IBM Decision Optimization on cloud. gap_best_obj = run_GAP_model(A, B, C) relaxed_best = run_GAP_model_with_Lagrangian_relaxation(A, B, C) # save the relaxed solution with get_environment().get_output_stream("solution.json") as fp: fp.write(json.dumps({"objectiveValue": relaxed_best}).encode('utf-8'))