Passing initial solution to quadratic and mixed integer problem

[Reinhold B]

Dear Prof. Löfberg,

as you clearly state, passing initial guesses for linear, quadratic and conic programming is typically neither necessary nor helpful.

But for integer programming it may help, but maybe less than I think.

Now I have a model consisting of a few quadratic constraints and some binary/integer constraints to obtain percentils of the previous squared variable.

You may find the simplified running example below.

As this combines a lot of integer logic I would like to pass an initial solution to the solver.

But as yalmip appends a normalized socp, the starting vector is not passed.

What would be necessary to provide a reasonable starting vector.

I already tried to turn off `convertconvexquad`, but then the quadratic constraint is not handled correctly.





Any comments on the problem are appreciated.


Are there different ways of expressing the problem, I could try?


In the end it is meant to optimize the input like a model predictive control based on the output of z


Best regards

Reinhold Bertram


example.m:

N = 20;
input = randn(1, N);

u = sdpvar(ones(1,N),ones(1,N));
y = sdpvar(ones(1,N),ones(1,N));

obj = 0;

ct = [ [u{:}] == input ];
ct = [ct, [y{:}] >= [u{:}].^2 ] ;

obj = obj + sum([y{:}]);

N_BINS = 10;

max_y = max([y{:}]);
min_y = min([y{:}]);
bin_width = (max_y - min_y)./N_BINS;

magnitude = sdpvar(N_BINS, 1);
edges = magnitude - bin_width/2;
edges(end+1,1) = magnitude(end) + bin_width/2;

MAX_MAGNITUDE = 100;
is_in_bin = binvar( repmat(N_BINS, 1, N), ones( 1, N));
bin_cnt = intvar(N_BINS, 1);

ct = [ct, [edges(1) == min_y]:['Edges 1'] ];
for n = 2:N_BINS
    ct = [ct, [edges(n) == (n-1) * bin_width + min_y]:['Edges ' num2str(n) ] ];
end
ct = [ct, [edges(N_BINS+1) == max_y]:['Edges ' num2str(N_BINS+1)] ];
          
ct = [ct, [0 <= magnitude <= MAX_MAGNITUDE]:['Bound magnitude' ] ];


for k = 1:N
    ct = [ct, [sum(is_in_bin{k}) == 1]:['sum isInBin ' num2str(k)] ];
    for n = 1:N_BINS
        ct = [ct,
            [ 0 <= y{k} <= MAX_MAGNITUDE]:['Bound Input']
            [implies(is_in_bin{k}(n), [edges(n) <= y{k} <= edges(n+1) ])]:['is in bin ' num2str(n) ':' num2str(k)];
            ];
    end
end

for n = 1:N_BINS
    sumInBin = 0;
    for k = 1:N
        sumInBin = sumInBin + is_in_bin{k}(n);
    end
    ct = [ct, [bin_cnt(n) == sumInBin]:['bin count'] ];
   
end

cum_probability = 100 * (1 - cumsum(bin_cnt) / N);

ct = [ct, [0 <= cum_probability <= 100]:['cum probability'] ];

CPF_MAP = [
    9
    17
    30
    ];
nMap = length(CPF_MAP);
z = sdpvar(nMap,1);

ct = [ct, [0 <= z <= 100]:['Bound Output']  ];

for ii = 1:nMap
    [~, loc] = min(abs(cum_probability - CPF_MAP(ii)));
    ct = [ct, [z(ii) ==  magnitude(loc)]:['Assign magnitude ' num2str(ii)] ];
end

%% Calculate initial guess

exp_y = input.^2;
[exp_bin_cnt, exp_magnitude] = hist(exp_y, N_BINS);
exp_bin_cnt = exp_bin_cnt(:);
exp_magnitude = exp_magnitude(:);

exp_bin_width = exp_magnitude(2) - exp_magnitude(1);
lb = exp_magnitude - exp_bin_width/2;
ub = exp_magnitude + exp_bin_width/2;
lb(1) = min(lb(1), min(exp_y));
ub(end) = max(ub(end), max(exp_y)) + eps;
mat_exp_y = repmat(exp_y, N_BINS, 1);
exp_is_in_bin = (lb <= mat_exp_y) & (mat_exp_y <= ub);

exp_z = nan(nMap, 1);
exp_cum_probability = 100 * (1 - cumsum(exp_bin_cnt) / N);
for ii = 1:nMap
    [~, idx] = min(abs(exp_cum_probability - CPF_MAP(ii)));
    exp_z(ii) = exp_magnitude(idx);
end

assign([u{:}], input);
assign([y{:}], exp_y);
assign(z, exp_z);
assign(magnitude, exp_magnitude);
assign([is_in_bin{:}], exp_is_in_bin);
assign(bin_cnt, exp_bin_cnt);

check(ct);

%%
options = sdpsettings('verbose', 2 ...
    , 'showprogress', 1 ...
    , 'solver', 'gurobi' ...
    , 'usex0', 1 ...
    , 'convertconvexquad', 1 ...
    );

diagnostic = optimize(ct, obj, options);

if diagnostic.problem == 0
    tc = matlab.unittest.TestCase.forInteractiveUse;
    tc.verifyEqual(value([u{:}]), input, 'AbsTol', 1e-6, 'u');
    tc.verifyEqual(value([y{:}]), exp_y, 'AbsTol', 1e-5, 'y');
    tc.verifyEqual(value(z), exp_z, 'AbsTol', 1e-6, 'z');
    tc.verifyEqual(value(magnitude), exp_magnitude, 'AbsTol', 1e-6, 'magnitude');
    tc.verifyEqual(value([is_in_bin{:}]), double(exp_is_in_bin), 'AbsTol', 1e-6, 'is_in_bin');
    tc.verifyEqual(value(bin_cnt), exp_bin_cnt, 'AbsTol', 1e-6, 'bin_cnt');
else
    disp(diagnostic);
end

[Johan Löfberg]

You have no way of passing initial values of variables introduced by operators such as implies. Some solvers accept partial initial guesses (supportsinitialNAN in definesolvers)

[Reinhold B]

As far as I understand append_normalized_socp, which is called during this example, it removes all given values for a possible x0.

And in the following even a partial initial guess is not possible, due to the missing initial values.

Would it be possible to append in append_normalized_socp, the x0 vector by NaNs to make the partial guess possible?

[Johan Löfberg]

I'll look into that

[Reinhold B]

Modifying append_normalized_socp.m with the attached patch, lets gurobi pick the MIPstart.

No idea of this works in all cases.

In the example case it speeds up the solution.

As you suggested in the blog entries mentioned above it does not help too much for my original problem.

Any suggestion how to reformulate the problem?

[Johan Löfberg]

I think a better way would be

if ~isempty(0)

        x0 = [x0;Ftemp*[1;x0]];

end

after Ftemp is defined. That way, all known initials are propagated correctly (assuming there are no nans, as those will kill the multiplication as 0*nan is nan. (A complete approach would be to first remove all zero columns of Ftemp corresponding to nan elements in x0)