# Issue in solving a large scale MIQP problem

I am solving a large scale MIQP optimisation problem at each step of a model predictive control problem. The problem description is as below. \begin{align} \min_{u} \quad (x_{k}&-x_\text{ref})^{T}Q(x_{k}-x_\text{ref}) + (P_{k}-P_{\text{ref},k})^{T}R(P_{k}-P_{\text{ref},k}) \\ \text{s.t. } x_{k+1}&=Ax_{k}+Bu_{k} \\ P_{k} &= \mathbb{I}^{T}u_{k}\\ u_{k} &\in \{0.25,0.50,0.75,1.0\}\\ \underline{x}&\leq x_{k}\leq \overline{x} \end{align} where $$Q$$ and $$R$$ are positive definite, $$\underline{x}$$ and $$\overline{x}$$ represent the lower and upper bounds of $$x$$, $$x_{ref}$$ and $$P_{ref}$$ are the reference values at each step $$k$$. In addition to that, $$P_{k} = \sum_{i} u_{k}^{i}$$ which is the sum of all decision variables at time step $$k$$.

The issue that I have is, the decision variable $$u$$ is in $$\mathbb{R}^{1000}$$, i.e. the problem involves a large number of integer variables which are not even binary.

I tried solving this optimisation problem at each iteration with Gurobi but was unable to solve at all. Thereafter, I contacted the Gurobi support centre and based on their suggestions, tried tweaking parameters (MIPGap, MIPHeuristics) to find at least a feasible solution. But it also did not improve the performance of the task. On the other hand, I am in need of finding a feasible solution within 60 seconds to match the real-world application.

I would really appreciate if someone could help me in the following problems.

1. Is this problem NP-hard? if so is there any way to solve this kind of a large problem based on the formulation shown above.
2. I am familiar with McCormick envelopes and big-M relaxations, but since the decision variables are not binary, is there any way I can apply those techniques here, I mean constraint relaxations?
3. Is there any powerful solver which I can try other than Gurobi?

Thank you.

• NP-hard. A problem with 1000 integer variables ain't chopped liver.In the worst case, it will never solve. Do you know a way to generate a feasible solution to provide to the solver? – Mark L. Stone Oct 5 at 1:54
• @MarkL.Stone I am not exactly aware of a feasible solution, somehow I need to find it by means of trial and error. I guess you are suggesting me to feed a feasible solution to the solver at the beginning, am I right ? If not could you please clarify your question a little bit ? – gayan_lanke Oct 5 at 2:12
• Since this is an MPC problem, what is the time horizon chosen? If you set $N$ to very low, then the problem should be fairly small. You'll get an aggressive control variable value, however it will be feasible. – Richard Oct 5 at 6:55
• Also, this is the post in the community forum of Gurobi related to this. Gurobi is able to find a solution, however it is slow at improving and it still has a 21.1% MIP gap after 2164 seconds. – Richard Oct 5 at 6:56
• @Richard, I tried with a prediction horizon $N$ of 3, still, cannot solve the problem. At the same time, it was me who posted on the Gurobi forum as well. However, the suggestions over there did not result in an attractive improvement. – gayan_lanke Oct 5 at 7:06

TL;DR: Some optimization problems are tough, and it requires a lot of work to get them solved.

1. Yes, this is NP-hard, but that does not say anything about whether or not it is easy to solve. Most MIPs are NP-hard, yet they are solved extremely frequently. My favorite treatment of this comes courtesy of Paul Rubin (see here).

2. You can always convert your integers into binaries. So if you have $$u \in \{0,1,2,3\}$$, you can convert this into $$y_k\in \{0,1\}$$, where $$u = \sum \limits_k ky_k$$, $$\sum \limits_k y_k = 1$$ and $$k\in\{0,1,2,3\}$$. Then you can apply the traditional McCormick relaxations and big-M formulations from the textbooks. Note that you can do McCormick relaxations also for a bound integer variable.

3. I work for Gurobi, so my obvious answer is "No" :) However, this being an MPC problem gives you actually quite a bit of structure. First off, I would try to work with shorter horizon lengths and see how that degrades performance. Next, I would look into relaxing $$u_k$$ and seeing what type of results you get. This is going to be the LP relaxation of the problem, and from the log file you posted on Gurobi's community forum it seems that the lower bound does not move much. So this may be a good place to start.

In general, you may have to create some specialized approaches. Note though that what you are really doing here is performance tuning, and following the comments of Yair Altman, you should always have a quantitative goal when performance tuning. So e.g. Gurobi does provide you a heuristic solution extremely quickly, however the MIP gap is fairly high. So you have to see which MIP gap would be acceptable for you and then add layers of algorithms on top to achieve this performance.

EDIT: thanks to Paul Rubin and Rob Pratt for spotting mistakes in my equations.

• Thank you for your response. I will take one by one and give my opinion. First of all, I started with the relaxed form of the QP problem where $u$ is not restricted to integers. It was very quick even for the case where $u$ $\in$ $\mathbb{R}^{1000}$ and the simulations results matched with my expectation for a tracking MPC scenario. Based on that, I realised that there is no flaw in the problem formulation. Thereafter, I tried the MIQP version where $u$ is only permitted to take integer values, and I could not solve the problem. – gayan_lanke Oct 5 at 7:53
• Considering 2), although that's a good idea, is there any benefit even if I reformulate to mimic a binary structure with $y_{k}$. Could you please clarify? Considering 3), I started with a large value for $N$ at the very beginning and then studied the performance and finally decided to choose $N=3$ to reduce the curse of dimensionality to a certain extent, and even for this, sometimes I might need to add a terminal cost or constraints to avoid stability issues. – gayan_lanke Oct 5 at 8:05
• Regarding > it seems that the lower bound does not move much. So this may be a good place to start. , actually, I did not get the exact idea. Would be happy if you could clarify it a bit. Finally, Ithe execution time is critical in my application, so I can only limit it for 60 seconds, and I even tried IterationLimit = 60 in Gurobi, but sometimes the solving process terminated at MIPGap around 70-80% which I cannot accept at all. – gayan_lanke Oct 5 at 8:05
• Thanks for the shout-out. I think there is a typo in (2), circa "where $y_k=k$", but I didn't want to edit it because I wasn't quite sure if you meant $u=\sum_k k y_k$ or something else. – prubin Oct 5 at 20:48
• For (2), you also need $\sum_k y_k =1$. – RobPratt Oct 6 at 12:32