I did not look too closely at your solution, but it seems that your original model was a convex QP problem. You made it a non-convex MINLP problem. I think you can find a solution that results in an BQP. First you need a set of binaries that tell you if ill>=2 (gama2) (via “classical” bigM methods), you also need a second set of binaries (b) that as soon as one turns 1 all following ones are also 1 (b(c,t+1)>=b(c,t)). For the second binary obviously needs to hold b(c,t) >= gama2(c,t). Moreover, b(c,t) can be used (with a bigM) to make qc (it seems you can also make ac) 0: ac(c,t) <= bigM*(1-b(c,t)). Your logic for the first condition is “ill(c,tp) =g= 2 * gama2(c,tp);”. That is questionable. ill(c,t) can be larger two and gama2(c,t) can still be 0. One does not always need to implement both direction of the equivalence some logic programming, the optimization may already do one direction automatically. I can’t easily figure this for your model.

Good luck,

-Michael