Non linear optimization with one constraint as exponential function

[Tam Karki]

Hi, 

I am looking to solve a non-linear optimization problem using YALMIP. For my problem, the MILNP solver should do the job, but not sure which one to use. 

The main concern is, to calculate the value of z, I have about 500 different exponential functions.

My objective is to Max: z

Subject to : 

z = exp((x-2)^2 + (y-5)^2) + exp((x-24)^2 + (y-10)^2) +......n , n is about 500 and all the terms are exponential functions

0<= x<= upper  bound

0 <= y <= upper bound

x and y both are non-negative.

I think this problem is non-convex. Not sure if YALMIP will handle it or not. Can anyone please help with where to start looking into to solve this problem optimally (global optimization would be the best)? 

Thanks,

Tam

[Johan Löfberg]

Somewhat smaller versions are solved rather easily by YALMIPs global solver (reaonable as the branching is performed in 2d so really pretty simple). Comparing that global solution with the solution from a local solver shows that the local solver typically is fairly ok, but often off by 10% or so

sdpvar x y

n = 50;

data = rand(n,2);

obj = -sum(exp(sum((repmat([x y],n,1)-data).^2,2)))

optimize(0 <= [x y] <= 1,obj)

obj

optimize(0 <= [x y] <= 1,obj,sdpsettings('solver','bmibnb'))

obj

[Tam Karki]

Thank you for your response. It helped a lot. I have a follow-up question on the YALMIP's global solver. 

Q1. Do you think YALMIP's global solver can handle the problem with a higher dimension (around 6d)?

Q2. What is the limitation on the size of the constraint when we use BMIBNB as a global solver (or in general YALMIP's global solvers)? Like I mentioned, I have constraints with exponential functions. For instance, z = exp((x-2)^2 + (y-5)^2) + 

exp((x-24)^2 + (y-10)^2) + 

exp((x-6)^2 + (y-2)^2) + 

exp((x-15)^2 + (y-5)^2) + 

exp((x-20)^2 + (y-5)^2) + 

exp((x-4)^2 + (y-10)^2) + 

exp((x-14)^2 + (y-5)^2) + 

exp((x-23)^2 + (y-10)^2) +

............+  exp((x-12)^2 + (y-5)^2) + 

exp((x-13)^2 + (y-10)^2) (there is a total of around 500-600 terms with exponential funtion in this constraints)

Thanks again!

Please let me know if I am not clear with my question.

Sincerely,

Tam

[Johan Löfberg]

Looks fairly easy

x = sdpvar(6,1);

z = randn(6,100);

objective = sum(exp(sum((repmat(x,1,100)-z).^2)));

optimize(-1 <= x <= 1,objective,sdpsettings('solver','bmibnb'))

* Starting YALMIP global branch & bound.

* Upper solver     : fmincon

* Lower solver     : MOSEK

* LP solver        : MOSEK

* -Extracting bounds from model

* -Perfoming root-node bound propagation

* -Feasible solution found by heuristics

* -Calling upper solver (found a solution!)

* -Branch-variables : 106

* -More root-node bound-propagation

* -Performing LP-based bound-propagation 

* -And some more root-node bound-propagation

* Starting the b&b process

 Node       Upper       Gap(%)       Lower     Open   Time

    1 :   9.32186E+06     0.00    9.32186E+06    0    10s  Terminated in bound propagation  

* Finished.  Cost: 9321856.1338 (lower bound: 9321856.1338, relative gap 0%)

* Termination with all nodes pruned 

* Timing: 11% spent in upper solver (1 problems solved)

*         0% spent in lower solver (0 problems solved)

*         81% spent in LP-based domain reduction (426 problems solved)

*         8% spent in upper heuristics (212 candidates tried)

ans = 

  struct with fields:

    yalmipversion: '20200930'

    matlabversion: '9.9.0.1524771 (R2020b) Update 2'

       yalmiptime: 1.0495

       solvertime: 10.2475

             info: 'Successfully solved (BMIBNB)'

          problem: 0

[Johan Löfberg]

but the problem does not have to be solved using bmibnb, since you easily can write it as a convex exponential cone problem (that's basically why it is solved immediately using bmibnb)

tisdag 23 mars 2021 kl. 16:19:59 UTC+1 skrev uttam...@gmail.com:

[Tam Karki]

Yeah, that totally makes sense. Thanks again! 

And what if I increase the bounds on my decision variables. x and y are actually the coordinates of a room. x and y both ranges from 0 to 1000 ft. Do you think the solver might increase the solving time significantly?

Sincerely,

Tam

[Johan Löfberg]

If you send it to bmibnb it might run into massive numerical issues with large bounds, as those exponentials essentially will be inf in the bounding which might cause issues. However, if it simply is written as the exponential cone program that it is, those bounds should not influence much.

Sounds like you are using the wrong scale if you expect stuff like exp((x-10)^2 with x is the range of thousands (since that exp will explode then)

[Tam Karki]

It was hard to write the original problem, so I wrote a demo one without thinking about scaling. Yeah, you are absolutely right, the constraint needs to be scaled. 

Thank you for your valuable comments. Will dig deeper into the exponential cone and BMIBNB solver. 

Sincerely,

Tam

[Johan Löfberg]

the expcone model is simply

t = sdpvar(1,100);

Model = [sum((repmat(x,1,100)-z).^2) <= t]

objective = sum(exp(t));

optimize(Model,objective)

which is solved using mosek if you have it,  otherwise it is just solved like any other nonlinear program. since it is convex they should have no issues with the model

[Tam Karki]

Sweet, thanks Prof. Löfberg!

Tam