IBMDecisionOptimization / docplex-examples

These samples demonstrate how to use the DOcplex library to model and solve optimization problems.
https://ibmdecisionoptimization.github.io/
Apache License 2.0
396 stars 229 forks source link

How to get daul_values #3

Closed keyfri closed 6 years ago

keyfri commented 6 years ago

When i tell Model.dual_values(cts), the output just print '0'.

` class Solution: def init(self, args):

general method parameter

    self.g_l = args.the_graph_length
    self.g_w = args.the_graph_width
    self.g_np = args.the_graph_length * args.the_graph_width
    self.g_toa = args.time_of_all
    self.g_max_u = self.g_np * self.g_toa
    self.g_m = args.num_UAV
    # abstract method parameter
    if args.the_graph_length % args.sparsity != 0 or args.the_graph_width % args.sparsity != 0 or args.time_of_all %  args.sparsity != 0:
        print('This game cannot sparsed!')
    else:
        self.ab_l = args.the_graph_length // args.sparsity
        self.ab_w = args.the_graph_width // args.sparsity
        self.ab_toa = args.time_of_all // args.sparsity
        self.ab_np = self.ab_l * self.ab_w
    # common parameter
    self.num_attack_zero = args.num_attack_zero
    self.a = args.arson_time_step
    self.d = args.possibility_of_same_node
    self.e = args.possibility_of_neighby_node

def clp(self, sub_set_of_attacker_strategies):

    #try:

        # Create model
        m = Model(name='CLP_model')
        # Variable
        flow = m.continuous_var_cube(self.g_np, self.g_np, self.g_toa + 1,
                           0.0, m.infinity, 'flow')
        mixPatrolVector = m.continuous_var_matrix(self.g_np, self.g_toa + 1,
                             0.0, float(self.g_m), 'mixPatrolVector')
        # Objective
        utility = m.continuous_var(- m.infinity, 0.0, 'utility')
        m.maximize(utility)

        # Constraint
        # Optimality constraint
        #optimality = m.add_range()
        m.optimality = []
        for index, path in enumerate(sub_set_of_attacker_strategies[0]):
            expr = m.linear_expr()
            uf, center_point = simple_utils(path, self.g_l, self.g_w,
                          self.num_attack_zero, self.a)
            expr.add_term(utility, -1.0)
            t = sub_set_of_attacker_strategies[1][index]
            for i in path:
                if i in center_point:
                    d = uf *self.d
                    expr.add_term(mixPatrolVector[(i, t)], d)
                    _, neighbor_point_arr = get_neighbor_point(get_axis_of_nodes(i,self.g_l), 1, self.g_l, self.g_w)
                    for neighbor_point in neighbor_point_arr:
                        for j in neighbor_point:
                            e = uf * self.e
                            expr.add_term(mixPatrolVector[(j, t)], e)
                t += 1
            #optimality.append(m.add(m.ge_constraint(expr, uf)))
            #optimality
            m.optimality.append(m.ge_constraint(expr, uf))
            m.add(m.ge_constraint(expr, uf))

        # out the flow
        for i in range(self.g_np):
            for t in range(1, self.g_toa):
                out_flow = m.linear_expr()
                out_flow.add_term(mixPatrolVector[(i, t)], -1)
                #out_flow.add_term(flow[(i, i, t)], 1)     # next step stay in  last point
                _, neighbor_point_arr = get_neighbor_point(get_axis_of_nodes(i, self.g_l), 1, self.g_l, self.g_w)
                for neighbor_point in neighbor_point_arr:
                    for j in neighbor_point:
                        out_flow.add_term(flow[(i, j, t)], 1)
                m.add(m.eq_constraint(out_flow, 0.0))

        # in the flow
        for i in range(self.g_np):
            for t in range(2, self.g_toa +1):
                in_flow = m.linear_expr()
                in_flow.add_term(mixPatrolVector[(i, t)], -1)
                #in_flow.add_term(flow[(i, i, t-1)], 1)    # last step stay in  next point
                _, neighbor_point_arr = get_neighbor_point(get_axis_of_nodes(i, self.g_l), 1, self.g_l, self.g_w)
                for neighbor_point in neighbor_point_arr:
                    for j in neighbor_point:
                        in_flow.add_term(flow[(j, i, t-1)], 1)
                m.add(m.eq_constraint(in_flow, 0.0))

        # flow
        for t in range(1, self.g_toa):
            start_flow = m.linear_expr()
            for i in range(self.g_np):
                start_flow.add_term(mixPatrolVector[(i, t)], 1)
            m.add(m.eq_constraint(start_flow, self.g_m))

        # Solve
        msol = m.solve()
        solution_value = m.objective_value
        flowValue = []
        for i in range(self.g_np):
            flow_j = []
            for j in range(self.g_np):
                flow_t = []
                for t in range(self.g_toa + 1):
                    flow_t.append(msol[flow[(i, j, t)]])
                flow_j.append(flow_t)
            flowValue.append(flow_j)

        mixPatrolValue = []
        for i in range(self.g_np):
            mixPatrol_t = []
            for t in range(self.g_toa + 1):
                mixPatrol_t.append(msol[mixPatrolVector[(i, t)]])
            mixPatrolValue.append(mixPatrol_t)
        #m.dual_value()
        #m.print_solution()
        duals  = m.dual_values(m.optimality)
        sub_set_of_attacker_strategies[2] = duals
        attacker_mix_patrol_strategy = sub_set_of_attacker_strategies
        m.end()
        return solution_value, flowValue, mixPatrolValue, attacker_mix_patrol_strategy

`

PhilippeCouronne commented 6 years ago

Hi I may have found an issue by looking into your code. When you build the "optimality" array, you add constraints that have not been added to the model. In this case, Docplex actually returns 0 without raising an exception (I think it should). The workaround is to store exactly those constraints that are added to the model, that is:

        ct = m.add(expr >= uf) # returns the added constraint
        m.optimality.append(ct) # store the added constraints.

In the code you sent, the constraints were duplicated; the call to dual values was receiving the copies unkwn to the model, hence the 0 value. Future versions of Docplex will make sure an exception is raised in this case.

However, as I could not run your code (extra functions, no data), this might not be the only issue. Let us know if the workaround works.

Regards.

Philippe Couronne (philippe.couronne@fr.ibm.com)

keyfri commented 6 years ago

Dear Philippe, Thank you for your reply. I have modified the code following your opinion. However, The dual_value still 0 value. So I share my program on GitHub https://github.com/keyfri/patrol_code/blob/master/src/game_solution.py . Looking forward to hearing from you!

sincerely, Jinpeng Han

PhilippeCouronne commented 6 years ago

Dear Jinpeng

Thanks for sharing your code, this was very helpful as I could try to run it. I encountered a runtime error at this line:

        #sub_set_of_attacker_strategies[2] = duals

where index 2 is out of bounds, so I commente dout this line. However, the computation of dual values is performed just above this line, so I was able to track this call to Cplex, and Cplex did return a list of zero values. As I recall, the dual value is the marginal impact on the cost of violating the constraint, but I could not relate the optimality constraints to the objective. You might need help from a mathematical expert in Linear Programming (I'm not) to investigate this further, sorry about that.

I had simplified your code to compute duals in:

duals = m.dual_values(m.optimality)

As Model.dual_values accepts sequences of constraints. I also wondered why you used factory functions (e.g. ge_constraint, le_constraint) rather than the overloaded Python operators such as >= , <=, ==

For example, instead of

            ct = m.add(m.ge_constraint(expr, uf))

You can write

            ct = m.ad(expr >= uf)

To summarize, with the modification, you code does call the Cplex dual values of the optimality constraint, and they are all zero.

PhilippeCouronne commented 6 years ago

Dear Jinpeng,

To investigate further your issue, I have modified your program to give names to optimality constraints, and exported the LP file. I then loaded this LP file into "cplex.exe", the interactive cplex optimizer, solved it and asked for the dualprices of all "optimality" constraints. Your issue is definitely a matehmatical one.

To add an optimality constraint with a name: ct = m.add(expr >= uf, "ct_optimality%d" % index) To export the model as a LP file: m.export_as_lp(basename="patrol") # relies on Python's tempfile module

Then uder cplex.exe, run the following commands

  1. read c:/temp/patrol.lp
  2. opt 3.display solution duals ct_optimality*

I get: Display dual values for which constraint(s): ct_optimality All dual prices matching 'ct_optimality' are 0.

Regards Philippe.