Very long solving time for MINLP

[Matthias Johansen]

Hi

I am working on a so called Security Constrained Unit Commitment problem with several generators in a 24 hour period. This means that the decision variable stating the power production is a nx24 matrix, where n is the number of generators. Same concept goes for the other decision variables, as you can see in the code below. To solve this problem I use Yalmip together with the Knitro solver. I am using the Knitro solver because this problem is a mixed integer non linear problem (quadratic?). The unit commitment problem becomes non linear when optimal power flow is added to the optimization model. 

There are 24 busses and 13 producing units.

My problem comes when I try to solve this model. It takes approximately 9000 seconds to solve this problem, and afterwards it comes up with a message saying "Optimal solution found".

Therefore the problem is not that the model is infeasible etc. The problem is the solving time.

Is this due to my Matlab code, or just a result of the non linearity coupled with binary decision variables?

In order to give some background for answers, the optimization code itself can be seen below.

In advance thank you for anwers. They are much appreciated.

%% Decision variables

% Creaing binary variables defined for turning on and off the engines

GenOnOff                   = binvar(NumGenUnits,h+1,'full'); % If 1 generation unit is on else off

% Creating variables that define the production of each unit

GenActivePowerProduction   = sdpvar(NumGenUnits,h,'full');

GenReactivePowerProduction = sdpvar(NumGenUnits,h,'full');

GenEfficiency              = sdpvar(NumGenUnits,h,'full');

% Creating binvar variables that defines each time a unit is started 

GenStart                   = binvar(NumGenUnits,h,'full');

GenStop                    = intvar(NumGenUnits,h,'full');

% Creating variables defining the SOC, Charge and discharge

SOC                        = sdpvar(NumBatUnits,h+1,'full');

Charge                     = sdpvar(NumBatUnits,h,'full');

DisCharge                  = sdpvar(NumBatUnits,h,'full');

BinCharge                  = binvar(NumBatUnits,h,'full');

% Creating Network values variables load flow constraints

CapBankON                  = binvar(NumCapBanks,h,'full');

BusVoltage                 = sdpvar(NumBus,h,'full');

VoltageAngles              = sdpvar(NumBus,h,'full');

RefBus                     = 1;%intvar(1,1);

%% Set up Constraint

constraints = [];

constraints = [constraints, GenOnOff(:,1) == 0];

    

for t = 1:24

    % Battery constraints

      constraints     = [constraints, ((1-BatDOD)*BatMaxCapacity) <= SOC(:,t+1) <= BatMaxCapacity];

      constraints     = [constraints, 0 <= Charge(:,t) <= (BatCRate*BatMaxCapacity.*BinCharge(:,t))./Sbase];

      constraints     = [constraints, 0 <= DisCharge(:,t) <= (BatCRate*BatMaxCapacity.*(1-BinCharge(:,t)))./Sbase];

      constraints     = [constraints, SOC(:,t+1) == SOC(:,t)-DisCharge(:,t).*Sbase+Charge(:,t).*Sbase];

    % Generator constraints

      constraints     = [constraints, abs(GenActivePowerProduction(:,t)+1i*GenReactivePowerProduction(:,t)) <= (GenMaxCapacity.*GenOnOff(:,t+1))./Sbase];

       constraints     = [constraints, (GenMinProduction.*GenOnOff(:,t+1))./Sbase <= GenActivePowerProduction(:,t)];

      constraints     = [constraints, GenOnOff(:,t+1)-GenOnOff(:,t) <= GenStart(:,t)];

      constraints     = [constraints, GenOnOff(:,t)-GenOnOff(:,t+1) == GenStop(:,t)];

      

       for unit = 1:NumGenUnits

           constraints     = [constraints, GenOnOff(unit,t+1:min(h,t+GenMinUpTime(unit))) >= GenStart(unit,t)];

           constraints     = [constraints, GenOnOff(unit,t+1:min(h,t+GenMinDownTime(unit))) <= 1-GenStop(unit,t)];

        end

      

    % Network constraints

      constraints = [constraints, (RenMatrix*RenProduction(t,:)')./Sbase+BatMatrix*DisCharge(:,t)-BatMatrix*Charge(:,t)+GenMatrix*(GenActivePowerProduction(:,t).*GenEfficiency(:,t))-PDemandMatrix(:,t) == real((BusVoltage(:,t).*exp(1i*VoltageAngles(:,t))).*conj(Ybus*                            (BusVoltage(:,t).*exp(1i*VoltageAngles(:,t)))))];

      constraints = [constraints, GenMatrix*GenReactivePowerProduction(:,t)+CapBankMatrix*(CapBank.*CapBankON(:,t)./Sbase)-QDemandMatrix(:,t) == imag((BusVoltage(:,t).*exp(1i*VoltageAngles(:,t))).*conj(Ybus*(BusVoltage(:,t).*exp(1i*VoltageAngles(:,t)))))];

      constraints = [constraints, 0.95 <= BusVoltage(:,t) <= 1.05];

      constraints = [constraints, VoltageAngles(RefBus,t) == 0];

      

    % N-1 criteria or 10% spinning reserve

     constraints = [constraints, max(sum((RenMatrix*RenProduction(t,:)')./Sbase+sum(BatMatrix*DisCharge(:,t)-BatMatrix*Charge(:,t))+sum(GenMatrix*(GenActivePowerProduction(:,t))))*0.1,max(GenActivePowerProduction(:,t))) <= sum((GenMaxCapacity./Sbase-          GenActivePowerProduction(:,t)).*GenOnOff(:,t+1))]; 

end

%% Define objective function

ObjGenProductionCosts = sum(sum(GenActivePowerProduction,2).*(GenMargProdCost+GenEmissionFactor));

ObjGenStartUpCosts    = sum(sum(GenStart,2).*GenStartUpCost);

Objective             = ObjGenProductionCosts + ObjGenStartUpCosts;

[Johan Löfberg]

without running code, it is impossible to decipher the model to see where the main complications are

It is not too uncommon that people destroy simple MILP-representable models using nonlinear constructs that can be expressed using linear combinatorial logic

[Matthias Johansen]

Hi Johan

Thanks for your answer. In order to run the code you need data from Excel files, and other scripts that contains the data used in the optimization script, so I think it would be too much for me to ask other people to run all this. If not please tell me.

Therefore the question is more general if MINLP (quadratic?) generally is very slow to solve. If it isn't, I will have to look at my code again.

The biggest complication, as far as I can see, is when the binary decision variable is used in the objective function, as well as the combination of non linearity and binary variables. I've tried to run the nonlinear part (AC power flow analysis) and the unit commitment model seperately. When doing that both models can be solved much faster. Therefore I suspect that the long solving time is due to the combination of nonlinearity and binary variables.

Even though i'm not nearly an expert in optimization models, I am 100% sure that a unit commitment with AC power flow is a nonlinear model.

Regards

Matthias

[Johan Löfberg]

Well, you seem to have exp in your model, so that sort of hints to me that you don't have a quadratic model.

Worst case with 24*13 binary variables, it could run for something like 10^67 times the age of the universe with a computer which can solve 1 billion relaxations per second ((2^(24*13))/(10^9*3600*24*365*14*10^9)). Absolutely impossible to say anything except that combinatorial / nonconvex nonlinear can have exponential complexity

Don't be shy to supply a self-contained package which runs easily

[Matthias Johansen]

Hi Johan

Okay, get your point :-)

Attached is a .zip file containing the used data and the Matlab script itself.

In row 9 and 10 the solver paths are set up. 

I don't know if you have access to Knitro. If you use another (maybe better) solver please let me know.

Again, thank you very much for your help.

Regards 

Matthias

[Johan Löfberg]

There are no low-hanging fruits here. This is a nasty nonlinear nonconvex model (polynomial/trigonometric) 

Also note that knitro is a local solver in the continuous space, hence the solutions you obtain are not necessarily optimal.

[Matthias Johansen]

So the long solving time is because of the nasty model. If a smaller solving time is extremely neccesary, I will need to remove some of the constraints that makes it this complex.

I know that Knitro is a local solver, but in my search after a proper solver I haven't found any better. The solver BARON is the one to use for a global optimal solution, but as far as I can read, it cannot manage my program, or maybe i'm wrong?.

 

[Johan Löfberg]

I don't see why it shouldn't work (i.e. run) with baron. I think it supports binary variables + polynomials and trigonometric function. It will perhaps run forever though

The long solution time is because your model is in theory impossible to solve to any reasonable degree of certified optimality. Practice may differ, but the outlook is not good.

[Matthias Johansen]

I understand that the model is too complex to find a good solution within a certain time frame. I will therefore try to reformulate the model into 2 different models. One MILP (unit commitment) and one NLP (AC power flow), where some of the values from the UC is used in the load flow analysis, and see if this helps.

Again thanks for your help.

[Mark L. Stone]

Per the latest BARON user manual (17.1.2, January 2, 2017)

In addition to multiplication and division, BARON can handle nonlinear functions that involve exp(x), ln(x), x for real α, and βx for real β. AIMMS/BARON, AMPL/BARON and
GAMS/BARON automatically handle |x| and xy, where x and y are variables; otherwise, suitable transformations discussed below can be used. Currently, there is no support for other functions, including the trigonometric functions sin(x), cos(x), etc.

[Johan Löfberg]

Ok. That's pretty strange, as it is trivial to implement those functions for a global branching procedure (easy enough for bmibnb in yalmip to support it)

[Mark L. Stone]

Here's what happens if you try using a trigonometric function with BARON under YALMIP

>> x = sdpvar;optimize(0<=x<=1,sin(x),sdpsettings('solver','baron','verbose',1,'debug',1))
Error using baron (line 179)
There was an error evaluating a user supplied function. Please examine the below error and correct your function.

Error using "callbaron/@(x)1*sin(x(1))" (line 0)
Undefined function 'sin' for input arguments of type 'barvec'.

Error in callbaron (line 79)
[x,fval,exitflag,info,allsol] = baron(obj,A,rl,ru,lb,ub,con,cl,cu,xtype,x0,opts);
Error in solvesdp (line 350)
    eval(['output = ' solver.call '(interfacedata);']);
Error in optimize (line 31)
[varargout{1:nargout}] = solvesdp(varargin{:});

[Mark L. Stone]

I think trig functions don't occur very often in the chemical engineering application which were the biggest application area as BARON was developed.