メインコンテンツへスキップ

Sensitivity analysis, quadratic constraints

回答済み
Benedicte Rimmereid
  • 編集済み
Conversationalist
First Question

I have a model that initially contained SOS constraints and quadratic constraints, on which I want to conduct a sensitivity analysis. First I used the fixed method to fix the SOS and then I tried to remove all the quadratic constraints after retrieving the information so I could query the linear constraints for constr.SARHSUp and constr.SARHSLow. See code below:

# Open the file in append mode, so existing content is preserved
with open("Results/sensitivity_analysisF.txt", "a") as file:

    # Remove quadratic constraints
    for qcon in fixed.getQConstrs():
        print(qcon)
        fixed.remove(qcon)
        fixed.update()

    print("hei")
    # Sensitivity analysis for non-quadratic constraints
    file.write("\nNon-Quadratic Constraint Sensitivity Analysis:\n")
    for constr in fixed.getConstrs():
        if fixed.NumQConstrs == 0:
            constr_name = constr.ConstrName
            constr_slack = constr.slack
            constr_dual_value = constr.Pi
            if constr_dual_value != 0:
                file.write(f"  Constraint: {constr_name}\n")
                file.write(f"  Slack: {constr_slack}\n")
                file.write(f"  Dual value (Shadow price): {constr_dual_value}\n")
                constr_rhs = constr.RHS
                file.write(f"  Right-hand side value: {constr_rhs}\n")
                constr_rhs_sensitivity_up = constr.SARHSUp
                constr_rhs_sensitivity_low = constr.SARHSLow
                file.write(f"  Right-hand side sensitivity (Increase): {constr_rhs_sensitivity_up}\n")
                file.write(f"  Right-hand side sensitivity (Decrease): {constr_rhs_sensitivity_low}\n")
Is there any known method I can use to retrieve the sensitivity analysis for the non-quadratic constraints?

関連するエンティティ

0

コメント

19件のコメント

日付 投票
  • Riley Clement
    Gurobi Staff Gurobi Staff

    Hi Benedicte,

    You re-optimized after removing the quadratic constraints right?

    - Riley

    0
  • Benedicte Rimmereid
    Conversationalist
    First Question

    Hi Riley, 

    Thank you for responding. The following code snippet might be more illustrative of my problem:

    file_path_sensitivity = os.path.join("Results", "SensitivityAdd.txt")
    fixed_q = read("fixed.lp")
    fixed_q.optimize()
    quad_constr = fixed_q.getQConstrs()
    if fixed_q.getConstrs() == None:
        print("Empty")
    else:
        print("Constraints to investigate")
        with open(file_path_sensitivity, "w") as file:
            file.write("\nConstraint sensitivity analysis:\n")
            for constr in fixed_q.getConstrs(): # for all non quadratic constraint
                if constr not in quad_constr:
                    constr_name = constr.ConstrName
                    print(constr_name)
                    constr_dual_value = constr.Pi
                    constr_slack = constr.Slack
                    if constr_dual_value != 0:
                        print(f"  Constraint: {constr_name}\n")
                        file.write(f"  Constraint: {constr_name}\n")
                        file.write(f"  Slack: {constr_slack}\n")
                        file.write(f"  Dual value (Shadow price): {constr_dual_value}\n")
                        constr_rhs = constr.RHS
                        file.write(f"  Right-hand side value: {constr_rhs}\n")
                        constr_rhs_sensitivity_up = constr.SARHSUp
                        constr_rhs_sensitivity_low = constr.SARHSLow
                        file.write(f"  Right-hand side sensitivity (Increase): {constr_rhs_sensitivity_up}\n")
                        file.write(f"  Right-hand side sensitivity (Decrease): {constr_rhs_sensitivity_low}\n")
     
    Here, I read a fixed version of my original problem and then found the quadratic constraints to avoid querying them for attributes only allowing for linear constraints. However, I get the following error:
    AttributeError: Unable to retrieve attribute 'Pi'
    0
  • Riley Clement
    Gurobi Staff Gurobi Staff

    Hi Benedicte,

    Can you try and set QCPdual parameter, as per Pi docs, and see how you go?

    - Riley

    0
  • Benedicte Rimmereid
    Conversationalist
    First Question

    Yes, I have tried that :/

     

     

    0
  • Riley Clement
    Gurobi Staff Gurobi Staff

    Maybe this example, based off qcp.py, will help?

    import gurobipy as gp
    from gurobipy import GRB

    # Create a new model
    m = gp.Model("qcp")

    # Create variables
    x = m.addVar(name="x", vtype="I")
    y = m.addVar(name="y")
    z = m.addVar(name="z")

    # Set objective: x
    obj = 1.0 * x
    m.setObjective(obj, GRB.MAXIMIZE)

    # Add constraint: x + y + z = 1
    c1 = m.addConstr(x + y + z == 1, "c0")

    # Add second-order cone: x^2 + y^2 <= z^2
    c2 = m.addConstr(x**2 + y**2 <= z**2, "qc0")

    # Add rotated cone: x^2 <= yz
    c3 = m.addConstr(x**2 <= y * z, "qc1")
    m.optimize()

    # end of qcp.py

    fixed = m.fixed()
    fixed.params.QCPDual=1
    fixed.optimize()

    for c in fixed.getConstrs():
        print(c.Pi)

    for qc in fixed.getQConstrs():
        print(qc.QCPi)

    Note the quadratic constraints need to be convex.

    - Riley

    0
  • Benedicte Rimmereid
    Conversationalist
    First Question

    Mhm, that results in:

    AttributeError: Unable to retrieve attribute 'Pi'

     

    0
  • Riley Clement
    Gurobi Staff Gurobi Staff
    • 編集済み

    What version of gurobipy are you using?

    You can verify with:

    print(gp.gurobi.version())
    0
  • Benedicte Rimmereid
    Conversationalist
    First Question

    (11,0,0) 

    :)

     

    0
  • Riley Clement
    Gurobi Staff Gurobi Staff

    Very strange, can you post the log output?

    0
  • Riley Clement
    Gurobi Staff Gurobi Staff

    Sorry, I meant the log output for the small example I posted.

    0
  • Benedicte Rimmereid
    Conversationalist
    First Question

    Ah okay :)


    Gurobi 11.0.0 (win64) logging started Tue Apr 16 14:54:47 2024

    Set parameter LogFile to value "gurobi_log_test.txt"
    Set parameter QCPDual to value 1
    Gurobi Optimizer version 11.0.0 build v11.0.0rc2 (win64 - Windows 11+.0 (22631.2))

    CPU model: AMD Ryzen 5 4500U with Radeon Graphics, instruction set [SSE2|AVX|AVX2]
    Thread count: 6 physical cores, 6 logical processors, using up to 6 threads

    Optimize a model with 1 rows, 3 columns and 3 nonzeros
    Model fingerprint: 0x31f5c2fc
    Model has 2 quadratic constraints
    Coefficient statistics:
      Matrix range     [1e+00, 1e+00]
      QMatrix range    [1e+00, 1e+00]
      Objective range  [1e+00, 1e+00]
      Bounds range     [0e+00, 0e+00]
      RHS range        [1e+00, 1e+00]
    Presolve removed 0 rows and 1 columns
    Presolve time: 0.03s
    Presolved: 2 rows, 2 columns, 3 nonzeros
    Presolved model has 1 second-order cone constraint
    Ordering time: 0.00s

    Barrier statistics:
     AA' NZ     : 1.000e+00
     Factor NZ  : 3.000e+00
     Factor Ops : 5.000e+00 (less than 1 second per iteration)
     Threads    : 1

                      Objective                Residual
    Iter       Primal          Dual         Primal    Dual     Compl     Time
       0  -0.00000000e+00 -0.00000000e+00  0.00e+00 1.00e-01  4.09e-02     0s
       1  -0.00000000e+00  4.99884378e-03  4.44e-16 3.69e-05  1.68e-03     0s
       2  -0.00000000e+00  9.33566787e-06  7.99e-15 4.06e-08  3.13e-06     0s
       3  -0.00000000e+00  1.74463971e-08  6.88e-15 4.47e-11  5.83e-09     0s

    Barrier solved model in 3 iterations and 0.06 seconds (0.00 work units)
    Optimal objective -0.00000000e+00

    Solving KKT system to obtain QCP duals...

    Optimize a model with 3 rows, 3 columns and 7 nonzeros
    Model fingerprint: 0x94eb7225
    Coefficient statistics:
      Matrix range     [4e-01, 1e+00]
      Objective range  [8e-01, 1e+00]
      QObjective range [2e+00, 2e+00]
      Bounds range     [0e+00, 0e+00]
      RHS range        [2e-01, 1e+00]
    Presolve removed 0 rows and 1 columns
    Presolve time: 0.08s
    Presolved: 3 rows, 2 columns, 6 nonzeros
    Ordering time: 0.00s

    Barrier statistics:
     AA' NZ     : 3.000e+00
     Factor NZ  : 6.000e+00
     Factor Ops : 1.400e+01 (less than 1 second per iteration)
     Threads    : 1

                      Objective                Residual
    Iter       Primal          Dual         Primal    Dual     Compl     Time
       0  -5.00000018e+05  5.00999975e+05  2.00e+03 0.00e+00  1.00e+06     0s
       1  -2.89435999e+03  5.45413893e+03  1.08e+02 0.00e+00  5.50e+04     0s
       2  -2.74230852e-01  2.64696428e+03  1.08e-04 0.00e+00  6.62e+02     0s
       3  -2.74155663e-01  3.80585952e+00  5.84e-08 0.00e+00  1.02e+00     0s
       4  -1.65218000e-01  6.40383121e-01  6.75e-09 0.00e+00  2.01e-01     0s
       5  -9.51900319e-05  4.55476129e-01  6.66e-15 0.00e+00  1.14e-01     0s
       6  -6.71885451e-06  4.50064159e-03  0.00e+00 5.55e-17  1.13e-03     0s
       7  -7.35822514e-12  4.70666831e-06  0.00e+00 0.00e+00  1.18e-06     0s
       8  -0.00000000e+00  4.71133366e-09  4.44e-16 1.11e-16  1.18e-09     0s

    Barrier solved model in 8 iterations and 0.08 seconds (0.00 work units)
    Optimal objective -0.00000000e+00

     

    0
  • Riley Clement
    Gurobi Staff Gurobi Staff

    Ok, that looks normal to me.

    I don't see why this approach would not work for your model unless your quadratic constraints were non-convex.  If you set NonConvex=1 on your fixed model and receive an error then your quadratic constraints are non-convex.

     

    0
  • Benedicte Rimmereid
    Conversationalist
    First Question

    Mhm, very wierd: 

    I tried this:

    fixed = model.fixed()
    fixed.params.QCPDual = 1 # TODO: back to 0
    quadratic_constraints = fixed.getQConstrs()
    fixed.params.NonConvex = 1
    fixed.update()    
    fixed.optimize()
    fixed.write("fixed.lp")
     
    And don't get any error, aka it should work? 
     
    0
  • Riley Clement
    Gurobi Staff Gurobi Staff

    That should work and write a QP.  If you want you can share the LP file here using dropbox, google drive etc

    0
  • Benedicte Rimmereid
    Conversationalist
    First Question

    https://drive.google.com/file/d/1UBCxGKIfO57T29zK4FCHM5cI4xnkB5nl/view?usp=sharing

    Thank you! 

    So to be clear, my end goal is to retrieve the SARHSUp/down-attribute

     

    0
  • Riley Clement
    Gurobi Staff Gurobi Staff

    Does this not work for you?

    m = gp.read("fixed.lp")
    m.params.QCPDual=1
    m.optimize()

    for c in m.getConstrs():
        print(c.Pi)

    for qc in m.getQConstrs():
        print(qc.QCPi)
    0
  • Benedicte Rimmereid
    Conversationalist
    First Question

    Yes, but not when I add SARHSUp :)

    0
  • Riley Clement
    Gurobi Staff Gurobi Staff

    Ok, so it sounds like this problem:

     AttributeError: Unable to retrieve attribute 'Pi' 

    is resolved.

    Looks like your original constraints and non PSD (so are not identified as convex) but Gurobi can translate them to convex SOCP constraints.  However these SOCP constraints prevent solutions from being basic, and SARHSUp and SARHSLow are only available for basic solutions.

    Perhaps there is a workaround where you remove all quadratic constraints and replace them with tangent planes at the optimal solution, but this would be quite a lot of work for you.

    - Riley

    0
  • Benedicte Rimmereid
    Conversationalist
    First Question

    Mhm I see,

    Thank you so much for the help.

    0

サインインしてコメントを残してください。

お探しのものを見つけられませんでしたか？

新規投稿