Infeasibility Issues / Optimal Control of Mixed Logical Dynamical Systems / MIQP

回答済み

Johannes Randenborgh

2023年03月02日 11:15

• 
• 
• 

Dear community,

first of all i'm new to gurobi with python. To study mixed integer quadratic programms (MIQP) with gurobi, i'm currently trying to implement an MIQP Problem based on Example 4.1 (System Dynamics) and Example 5.1 (Optimal Control) published by Bemporad, Morari 1999 Control of Systems Integrating logic, dynamics, and constraints, Automatica 35, p. 407-427.

The system dynamics, a piecewise linear time invaritant system, can be transformed into a MIQP with linear inequality constraints and quadratic, convex, cost. The MIQP looks as follows:

V ist the optimization variable ( addMVar(): a stacked vector of input (shape: nu,1), state (shape: nz,1) and binary (shape: nd,1) variables with GRB.CONTINUOUS and GRB.BINARY); S1, S2, S3, F1, F2, F3 are constant matrices, x0 is the starting state of the system, a constant.

Solving the optimization problem, gurobi informs that the problem is infeasible. From my point of view, infeasiblity may only result from erros in the inequality constraints. However, i was able to check the inequality constraints by sucessfully testing them with a possible V, x0, what would be a contradiction, right? Thus, i believe I have not properly implemented the optimization problem with gurobi.py - syntax errors?. You can see my code below.  The given equality constraint is not brought into an inequality constraint formulation to keep it simple. This constraint results from the logic of the system dynamics and is necessary.

try:
    # create model
    m = gp.Model(' ex41')

    # create variables
    V = m.addMVar(shape=(N*(nu + nd + nz), 1), vtype = GRB.CONTINUOUS, name='V')
    V[N*nu:N*(nu+nd)].VType = GRB.BINARY # define binary variables

    d = V[N*nu:N*(nd+nu)] # extract binary variables
    
    obj = V.T @ S1 @ V + 2 * (S2 + x0.T @ S3) @ V

    m.setObjective(obj, GRB.MINIMIZE)

    # define constraints
    # F1 V <= F2 + F3 x0
    m.addConstr( F1 @ V <= F2 + F3 @ x0, name='ineq')
 
    # set condition - only one set active 
    m.addConstr( np.kron(np.eye(N), np.ones((1,ns))) @ d == np.ones((N,1)), name='xor') # sum(di) = 1

    m.setParam('DualReductions', 0)
    m.setParam('IntFeasTol', 1e-5) # 1e-5 is default
    m.update()
    m.optimize()
    m.printStats()

    m.computeIIS()
    m.write('model.ilp')

    if m.status == GRB.INFEASIBLE:
        m.feasRelaxS(1, False, False, True)
        m.optimize()

except gp.GurobiError as e:
    print('Error code ' + str(e.errno) + ": " + str(e))

except AttributeError:
    print('Encountered an attribute error')

The solver output is given here:

Optimize a model with 58 rows, 14 columns and 136 nonzeros
Model fingerprint: 0x9fb4bfbb
Model has 16 quadratic objective terms
Variable types: 10 continuous, 4 integer (4 binary)
Coefficient statistics:
  Matrix range     [4e-01, 1e+03]
  Objective range  [2e-02, 2e-02]
  QObjective range [2e-02, 4e+00]
  Bounds range     [0e+00, 0e+00]
  RHS range        [1e+00, 1e+03]
Presolve removed 25 rows and 0 columns
Presolve time: 0.00s

Explored 0 nodes (0 simplex iterations) in 0.00 seconds (0.00 work units)
Thread count was 1 (of 8 available processors)

Solution count 0

Model is infeasible
Best objective -, best bound -, gap -

Statistics for model ex41:
  Linear constraint matrix    : 58 Constrs, 14 Vars, 136 NZs
  Variable types              : 10 Continuous,
4 Integer (4 Binary)
  Matrix coefficient range    : [ 0.4, 1014.66 ]
  Objective coefficient range : [ 0.02, 0.02 ]
  Variable bound range        : [ 0, 0 ]
  RHS coefficient range       : [ 1, 1020.12 ]

Computing the IIS model and looking at 'model.ils' do not properly help me to identify the error.

Relaxing constraints works. However, the results are far away of the constraints and cannot be used.

Would it be possible to check on my syntax and other possible errors? In case you need the code for the matrices (F1, F2, ...), please let me know.

Bests,

関連

• Python
• Infeasible

0

• 

  Jaromił Najman

  • Gurobi Staff

  2023年03月02日 13:21

  Hi Johannes,

  You could write the model you construct to a human-readable LP file via the write method.

  [...]
  m.update()
  
  m.write("myModel.lp")
  
  m.optimize()
  m.printStats()
  [...]

  You can open the file \(\texttt{myModel.lp}\) in any standard Text Editor and analyze whether the model looks as expected.

  One thing that comes to mind is that the default lower bound value for optimization variables in Gurobi is \(0\). If your variables are allowed to take negative values, then you have to explicitly state so when constructing the variables.

  V = m.addMVar(shape=(N*(nu + nd + nz), 1), lb=-GRB.INFINITY, vtype = GRB.CONTINUOUS, name='V')

  You could also share the content of the LP file here if further help is needed.

  Best regards, 
  Jaromił

   

  0
• 

  Johannes Randenborgh

  2023年03月02日 18:00

  • 
  • 
  • 

  Hi Jaromił,

  thanks for the immediate reply.

  **Solution**: I forgot so set lb= - GRB.INFINITY. The system dynamics force the state to alternate between positive and negative values.

  Bests,

  Johannes

  0

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