Optimization with precalculated arrays: How to get values from the array? If it is not possible, how can I approache the problem in differently?

[Konstantin Römer]

Dear all,

i got a problem regarding an optimization with yalmip. I try to minimize the cost of a product and use precalculated values for that. How do I approach the problem of indexing and false logical constraints?
I first tried to build a constraint directly out of the array: Constraints = [Constraints, NPV(PV, BAT) >= 0];
Yalmip throws here as expected the error: "Function 'subsindex' is not defined for values of class 'sdpvar'."

So my first question is, can i some how used the constraint as mentioned above or is it not possible?

Second, i tried to transform the full array into a set of contraints.
for n = 1:10
    for m = 1:10
        if NPV(n,m) < 0
            Constraints = [Constraints, implies(PV == n, BAT ~= m)];
        end
    end
end
In addition, I change the objective function into:
minCost = (PV * costPV + BAT * costBAT + PV * costINV);

However, it would be very nice if i could rewrite the expression as follows:
minCost = (PV * costPV + BAT * costBAT + INVsize(PV, BAT) * costINV);

Thank you for your support. I'm just a beginner in Matlab and not very expierienced.
Best

[Johan Löfberg]

MATLAB does not support overloading of indexing a double with user class

You can get by this by defining data = sdpvar(10,10,'full'), add the constraint data == NPV, and then reference the matrix data instead.

Typically though, this overloaded indexing in a general scenario is an expensive thing in terms of the resulting MILP model (and I haven't really tried to optimize it), and you might gain from doing this manually. I guess one model is

d=binvar(10,10);
thevalue_at_PV_BAT = sum(sum(NPV.*d))
Model = [sum(sum(d))==1,sum(sum(d*(1:10)'))==PV,sum(sum((1:10)*d))==BAT]

[Konstantin Römer]

Thank you for the fast response. I was able to implement the optimization.

I implemented the following
dimn = 10;
dimm = 10;
d=binvar(dimn, dimm);
npv = sum(sum(NPV.*d));

The resulting constraint is
Constraints = [Constraints, npv >= 0];

I concluded the following objective function:
minCost = (PV * costPV + BAT * costBAT + inv(PV,BAT) * costINV + objc(PV, BAT));

where objc(PV,BAT) is a fixed value at the point (PV, BAT)
Did i implement that correct?
Thank you for your support!

[Johan Löfberg]

No, this is something completely different. Since d is a completely unconstrained matrix, the constraint npv>=0 is entirely redundant as d=0 is feasible

In your text above, you only wanted NPV(PV,BAT)>=0, hence d and the constraints connecting the binary pattern d with the indicies PV and BAT, and the constraint that only one element in d is used, are crucial. Also, the flag 'full' is not something I added randomly to d...

and objc(PV,BAT), wasn't the whole problem the fact that you cannot index a constant matrix

[Konstantin Römer]

I'm Sorry, I totally forgot to add that i constrained d as suggested:
Constraints = [Constraints, sum(sum(d)) == 1, sum(sum(d*(1:dimn)'))==PV,sum(sum((1:dimm)*d))==BAT];

However, the simulation results always into the same value. I attached the full code.

[Johan Löfberg]

Code doesn’t make sense. You are mixing the manual modelling strategy of indexing with decision variables (i.e. using d) with the other approach based on rewriting to indexing in a decision variable

[Konstantin Römer]

Thanks for your response. I tried to change it but im not sure how you thought it would be better. I realy tried to implement as suggested but im not sure if it correct now.
One issue that occures is that i recieve the same values for the decision variables PV and BAT?

Thank you for your support!

Best
Konstantin

[Johan Löfberg]

One weird thing with your model is that you want to ensure NPV(PV,BAT)>=0 and do that by a constraint. Why do you even have that possibility in the data?

The indexing model is correct though, so if you get PV=BAT, it has to be something inherent in the rest of your model, or something specific to your data (have you tried random?)

This example problem shows that the idea works

dimn = 10;

dimm = 20;

intvar PV BAT

d = binvar(dimm,dimn);

data = randn(dimm,dimn)*100;

Model = [sum(sum(d)) == 1, sum(sum(d*(1:dimn)'))==PV ,sum(sum((1:dimm)*d))==BAT];

objective = sum(sum(data.*d)) + PV^2+BAT^2

optimize(Model,objective)

% Optimal objective from optimal indicies

data(value(BAT),value(PV))+value(PV)^2+value(BAT)^2

% Compare with brute-force

[pv,bat] = meshgrid(1:dimn,1:dimm);

Z = pv.^2+bat.^2+data;

min(min(Z))

[Johan Löfberg]

BTW, for the particular problem you posted now, there is absolutely no reason to use integer programming (or any solver at all),. You will never beat a trivial double for-loop which simply goes through all cases and computes the objective, followed by picking the best. You don't have any constraints in your model.