375 lines
13 KiB
Mathematica
375 lines
13 KiB
Mathematica
|
function [model,keptvariablesIndex] = eliminatevariables(model,removedparameters,value,parameters)
|
||
|
|
|
||
|
|
if isempty(removedparameters)
|
||
|
|
keptvariablesIndex = 1:length(model.c);
|
||
|
|
return
|
||
|
|
end
|
||
|
|
|
||
|
|
keptvariablesIndex = 1:length(model.c);
|
||
|
|
|
||
|
|
rmvmonoms = model.rmvmonoms;
|
||
|
|
newmonomtable = model.newmonomtable;
|
||
|
|
|
||
|
|
% Assign evaluation-based values
|
||
|
|
if length(model.evalParameters > 0)
|
||
|
|
alls = [removedparameters(:)' model.evalParameters(:)'];
|
||
|
|
removedP = [];
|
||
|
|
removedE = [];
|
||
|
|
for j = 1:length(removedparameters)
|
||
|
|
for i = 1:length(model.evalMap)
|
||
|
|
if isequal(model.evalMap{i}.variableIndex,removedparameters(j))
|
||
|
|
p = value(j);
|
||
|
|
index = find(model.evalMap{i}.computes == alls);
|
||
|
|
value(index) = feval(model.evalMap{i}.fcn,p);
|
||
|
|
removedP = [removedP model.evalMap{i}.computes];
|
||
|
|
removedE = [removedE i];
|
||
|
|
end
|
||
|
|
end
|
||
|
|
end
|
||
|
|
model.evalVariables(removedE)=[];
|
||
|
|
model.evalMap = {model.evalMap{setdiff(1:length(model.evalMap),removedE)}};
|
||
|
|
model.evalParameters = setdiff(model.evalParameters,removedP);
|
||
|
|
end
|
||
|
|
|
||
|
|
% Evaluate fixed monomial terms
|
||
|
|
aux = model.precalc.aux;
|
||
|
|
if ~isempty(model.precalc.jj1)
|
||
|
|
z = value(model.precalc.jj1).^rmvmonoms(model.precalc.index1);
|
||
|
|
aux(model.precalc.index1) = z;
|
||
|
|
end
|
||
|
|
% Assign simple linear terms
|
||
|
|
if ~isempty(model.precalc.jj2)
|
||
|
|
aux(model.precalc.index2) = value(model.precalc.jj2);
|
||
|
|
end
|
||
|
|
|
||
|
|
monomvalue = prod(aux,2);
|
||
|
|
|
||
|
|
removethese = model.removethese;
|
||
|
|
keepingthese = model.keepingthese;
|
||
|
|
|
||
|
|
% Look for constraints which evaluate to inf + a'*x + b >= 0
|
||
|
|
infinite_monoms = removethese(isinf(monomvalue(removethese)));
|
||
|
|
if ~isempty(infinite_monoms)
|
||
|
|
if model.K.l > 0
|
||
|
|
A = model.F_struc(model.K.f + 1:model.K.f + model.K.l,1 + infinite_monoms);
|
||
|
|
redundant = find(all(A,2));
|
||
|
|
model.F_struc(model.K.f + redundant,:)=[];
|
||
|
|
model.K.l = model.K.l - length(redundant);
|
||
|
|
end
|
||
|
|
% If that infinite variable is used nowhere else, simply zero it so we
|
||
|
|
% don't get any inf*0 expressions as these leads to nan
|
||
|
|
if ~any(model.F_struc(:,1+infinite_monoms))
|
||
|
|
monomvalue(infinite_monoms)=0;
|
||
|
|
end
|
||
|
|
end
|
||
|
|
|
||
|
|
value = monomvalue(removethese);
|
||
|
|
monomgain = monomvalue(keepingthese);
|
||
|
|
|
||
|
|
if ~isempty(model.F_struc)
|
||
|
|
model.F_struc(:,1) = model.F_struc(:,1)+model.F_struc(:,1+removethese)*value;
|
||
|
|
model.F_struc(:,1+removethese) = [];
|
||
|
|
model.F_struc = model.F_struc*diag(sparse([1;monomgain]));
|
||
|
|
end
|
||
|
|
infeasible = 0;
|
||
|
|
if model.K.f > 0
|
||
|
|
candidates = find(~any(model.F_struc(1:model.K.f,2:end),2));
|
||
|
|
%candidates = find(sum(abs(model.F_struc(1:model.K.f,2:end)),2) == 0);
|
||
|
|
if ~isempty(candidates)
|
||
|
|
% infeasibles = find(model.F_struc(candidates,1)~=0);
|
||
|
|
if find(model.F_struc(candidates,1)~=0,1)%;~isempty(infeasibles)
|
||
|
|
model.infeasible = 1;
|
||
|
|
return
|
||
|
|
else
|
||
|
|
model.F_struc(candidates,:) = [];
|
||
|
|
model.K.f = model.K.f - length(candidates);
|
||
|
|
end
|
||
|
|
end
|
||
|
|
end
|
||
|
|
if model.K.l > 0
|
||
|
|
% Find infeasible constraints
|
||
|
|
candidates = find(~any(model.F_struc(model.K.f + (1:model.K.l),2:end),2));
|
||
|
|
if ~isempty(candidates)
|
||
|
|
if any(model.F_struc(model.K.f + candidates,1)<-1e-14)
|
||
|
|
model.infeasible = 1;
|
||
|
|
return
|
||
|
|
else
|
||
|
|
model.F_struc(model.K.f + candidates,:) = [];
|
||
|
|
model.K.l = model.K.l - length(candidates);
|
||
|
|
end
|
||
|
|
end
|
||
|
|
% Find constraints f(x) <= inf and remove these
|
||
|
|
candidates = find(isinf(model.F_struc(model.K.f + (1:model.K.l),1)) & (model.F_struc(model.K.f + (1:model.K.l),1))>0);
|
||
|
|
if ~isempty(candidates)
|
||
|
|
model.F_struc(model.K.f + candidates,:) = [];
|
||
|
|
model.K.l = model.K.l - length(candidates);
|
||
|
|
end
|
||
|
|
end
|
||
|
|
if model.K.q(1) > 0
|
||
|
|
removeqs = [];
|
||
|
|
removeRows = [];
|
||
|
|
top = model.K.f + model.K.l + 1;
|
||
|
|
F_struc = model.F_struc(top:top+sum(model.K.q)-1,:);
|
||
|
|
top = 1;
|
||
|
|
if all(any(F_struc(:,2:end),2))
|
||
|
|
% There is still something on every row in all SOCPs, we don't have
|
||
|
|
% to search for silly SOCPS ||0|| <= constant
|
||
|
|
else
|
||
|
|
for i = 1:length(model.K.q)
|
||
|
|
rows = top:top+model.K.q(i)-1;
|
||
|
|
v = F_struc(rows,:);
|
||
|
|
if nnz(v(:,2:end))==0
|
||
|
|
if norm(v(2:end,1)) > v(1,1)
|
||
|
|
model.infeasible = 1;
|
||
|
|
return
|
||
|
|
else
|
||
|
|
removeqs = [removeqs;i];
|
||
|
|
removeRows = [removeRows;model.K.f+model.K.l+rows];
|
||
|
|
end
|
||
|
|
end
|
||
|
|
top = top + model.K.q(i);
|
||
|
|
end
|
||
|
|
model.K.q(removeqs)=[];
|
||
|
|
model.F_struc(removeRows,:)=[];
|
||
|
|
if isempty(model.K.q)
|
||
|
|
model.K.q = 0;
|
||
|
|
end
|
||
|
|
end
|
||
|
|
end
|
||
|
|
if model.K.s(1) > 0
|
||
|
|
% Nonlinear semidefinite program with parameter
|
||
|
|
top = model.K.f + model.K.l + sum(model.K.q) + 1;
|
||
|
|
removeqs = [];
|
||
|
|
removeRows = [];
|
||
|
|
for i = 1:length(model.K.s)
|
||
|
|
n = model.K.s(i);
|
||
|
|
rows = top:top+n^2-1;
|
||
|
|
v = model.F_struc(rows,:);
|
||
|
|
if nnz(v)==0
|
||
|
|
removeqs = [removeqs;i];
|
||
|
|
removeRows = [removeRows;rows];
|
||
|
|
elseif nnz(v(:,2:end))==0
|
||
|
|
Q = reshape(v(:,1),n,n);
|
||
|
|
used = find(any(Q));Qred=Q(:,used);Qred = Qred(used,:);
|
||
|
|
[~,p] = chol(Qred);
|
||
|
|
if p
|
||
|
|
model.infeasible = 1;
|
||
|
|
return
|
||
|
|
else
|
||
|
|
removeqs = [removeqs;i];
|
||
|
|
removeRows = [removeRows;rows];
|
||
|
|
end
|
||
|
|
end
|
||
|
|
top = top + n^2;
|
||
|
|
end
|
||
|
|
model.K.s(removeqs)=[];
|
||
|
|
model.F_struc(removeRows,:)=[];
|
||
|
|
if isempty(model.K.s)
|
||
|
|
model.K.s = 0;
|
||
|
|
end
|
||
|
|
end
|
||
|
|
|
||
|
|
model.f = model.f + model.c(removethese)'*value;
|
||
|
|
|
||
|
|
model.c(removethese)=[];
|
||
|
|
if nnz(model.Q)>0
|
||
|
|
model.c = model.c + 2*model.Q(keepingthese,removethese)*value;
|
||
|
|
model.f = model.f + value'*model.Q(removethese,removethese)*value;
|
||
|
|
end
|
||
|
|
|
||
|
|
if nnz(model.Q)==0
|
||
|
|
model.Q = spalloc(length(model.c),length(model.c),0);
|
||
|
|
else
|
||
|
|
model.Q(removethese,:) = [];
|
||
|
|
model.Q(:,removethese) = [];
|
||
|
|
end
|
||
|
|
model.c = model.c.*monomgain;
|
||
|
|
keptvariablesIndex(removethese) = [];
|
||
|
|
|
||
|
|
model.lb(removethese)=[];
|
||
|
|
model.ub(removethese)=[];
|
||
|
|
|
||
|
|
newmonomtable(:,removethese) = [];
|
||
|
|
newmonomtable(removethese,:) = [];
|
||
|
|
|
||
|
|
if ~isequal(newmonomtable,model.precalc.newmonomtable)%~isempty(removethese)
|
||
|
|
skipped = [];
|
||
|
|
alreadyAdded = zeros(1,size(newmonomtable,1));
|
||
|
|
[ii,jj,kk,skipped] = stableunique(newmonomtable*model.hashCache(1:size(newmonomtable,2)));
|
||
|
|
%[ii,jj,kk,skipped] = stableunique(newmonomtable*gen_rand_hash(0,size(newmonomtable,2),1));
|
||
|
|
S = sparse(kk,1:length(kk),1);
|
||
|
|
model.precalc.S = S;
|
||
|
|
model.precalc.skipped = skipped;
|
||
|
|
model.precalc.newmonomtable = newmonomtable;
|
||
|
|
model.precalc.blkOneS = blkdiag(1,S');
|
||
|
|
else
|
||
|
|
S = model.precalc.S;
|
||
|
|
skipped = model.precalc.skipped;
|
||
|
|
end
|
||
|
|
model.c = S*model.c;
|
||
|
|
if ~isempty(model.F_struc)
|
||
|
|
model.F_struc = model.F_struc*model.precalc.blkOneS;
|
||
|
|
end
|
||
|
|
|
||
|
|
model.lb(skipped) = [];
|
||
|
|
model.ub(skipped) = [];
|
||
|
|
newmonomtable(skipped,:) = [];
|
||
|
|
newmonomtable(:,skipped) = [];
|
||
|
|
|
||
|
|
|
||
|
|
if nnz(model.Q)==0
|
||
|
|
model.Q = spalloc(length(model.c),length(model.c),0);
|
||
|
|
else
|
||
|
|
model.Q(:,skipped)=[];
|
||
|
|
model.Q(skipped,:)=[];
|
||
|
|
end
|
||
|
|
|
||
|
|
keptvariablesIndex(skipped) = [];
|
||
|
|
|
||
|
|
model.monomtable = newmonomtable;
|
||
|
|
model = compressModel(model);
|
||
|
|
|
||
|
|
x0wasempty = isempty(model.x0);
|
||
|
|
model.x0 = zeros(length(model.c),1);
|
||
|
|
|
||
|
|
% Try to reduce to QP
|
||
|
|
[model,keptvariablesIndex,newmonomtable] = setupQuadratics(model,keptvariablesIndex,newmonomtable);
|
||
|
|
|
||
|
|
if nnz(model.Q) > 0
|
||
|
|
if model.solver.objective.quadratic.convex == 0 & model.solver.objective.quadratic.nonconvex == 0
|
||
|
|
error('The objective instantiates as a quadratic after fixing parameters, but this is not directly supported by the solver. YALMIP will not reformulate models if they structurally change in call to optimizer. A typical trick to circumvent this is to define a new set of variable e, use the quadratic function e''e, and add an equality constraint e = something. The SOCP formulation can then be done a-priori by YALMIP.');
|
||
|
|
end
|
||
|
|
end
|
||
|
|
|
||
|
|
if x0wasempty
|
||
|
|
model.x0 = [];
|
||
|
|
end
|
||
|
|
|
||
|
|
% Remap indicies
|
||
|
|
if ~isempty(model.integer_variables)
|
||
|
|
temp=ismember(keptvariablesIndex,model.integer_variables);
|
||
|
|
model.integer_variables = find(temp);
|
||
|
|
end
|
||
|
|
if ~isempty(model.binary_variables)
|
||
|
|
temp=ismember(keptvariablesIndex,model.binary_variables);
|
||
|
|
model.binary_variables = find(temp);
|
||
|
|
end
|
||
|
|
if ~isempty(model.semicont_variables)
|
||
|
|
temp=ismember(keptvariablesIndex,model.semicont_variables);
|
||
|
|
model.semicont_variables = find(temp);
|
||
|
|
end
|
||
|
|
if ~isempty(model.aux_variables)
|
||
|
|
temp=ismember(keptvariablesIndex,model.aux_variables);
|
||
|
|
model.aux_variables = find(temp);
|
||
|
|
end
|
||
|
|
if ~isempty(model.K.sos)
|
||
|
|
if ~isempty(model.K.sos.variables)
|
||
|
|
for i = 1:length(model.K.sos.variables)
|
||
|
|
temp=ismember(keptvariablesIndex,model.K.sos.variables{i});
|
||
|
|
model.K.sos.variables{i} = find(temp);
|
||
|
|
end
|
||
|
|
end
|
||
|
|
end
|
||
|
|
|
||
|
|
model.used_variables = model.used_variables(keptvariablesIndex);
|
||
|
|
|
||
|
|
% Check if there are remaining strange terms. This occurs in #152
|
||
|
|
% FIXME: Use code above recursively instead...
|
||
|
|
if ~model.solver.objective.sigmonial & any(model.variabletype == 4)
|
||
|
|
% Bugger. There are signomial terms left, despite elimination, and the
|
||
|
|
% solver does not handle this. YALMIP has introduced an intermediate
|
||
|
|
% variable which is a nasty function of the parameter.
|
||
|
|
signomials = find(model.variabletype == 4);
|
||
|
|
involved = [];
|
||
|
|
for i = 1:length(signomials)
|
||
|
|
m = model.monomtable(signomials(i),:);
|
||
|
|
involved = [involved;find(m ~= fix(m) | m < 0)];
|
||
|
|
end
|
||
|
|
involved = unique(involved);
|
||
|
|
[lb,ub] = findulb(model.F_struc,model.K,model.lb,model.ub);
|
||
|
|
if all(lb(involved) == ub(involved))
|
||
|
|
% Now add equality constraints to enforce
|
||
|
|
for i = signomials
|
||
|
|
m = model.monomtable(i,:);
|
||
|
|
involved = find(m ~= fix(m) | m < 0);
|
||
|
|
gain = lb(involved).^m(involved);
|
||
|
|
s = zeros(1,size(model.F_struc,2));
|
||
|
|
multiplies = setdiff(find(m),involved);
|
||
|
|
s(i+1) = 1;
|
||
|
|
s(multiplies+1) = -gain;
|
||
|
|
model.F_struc = [s;model.F_struc];
|
||
|
|
model.K.f = model.K.f + 1;
|
||
|
|
end
|
||
|
|
else
|
||
|
|
error('Did not manage to instatiate model. Complicating terms remaining');
|
||
|
|
end
|
||
|
|
end
|
||
|
|
|
||
|
|
function model = compressModel(model)
|
||
|
|
model.variabletype = zeros(size(model.monomtable,1),1)';
|
||
|
|
nonlinear = sum(model.monomtable,2)~=1 | sum(model.monomtable~=0,2)~=1;
|
||
|
|
if ~isempty(nonlinear)
|
||
|
|
model.variabletype(nonlinear) = 3;
|
||
|
|
quadratic = sum(model.monomtable,2)==2;
|
||
|
|
model.variabletype(quadratic) = 2;
|
||
|
|
bilinear = max(model.monomtable,[],2)<=1;
|
||
|
|
model.variabletype(bilinear & quadratic) = 1;
|
||
|
|
sigmonial = any(0>model.monomtable,2) | any(model.monomtable-fix(model.monomtable),2);
|
||
|
|
model.variabletype(sigmonial) = 4;
|
||
|
|
end
|
||
|
|
|
||
|
|
|
||
|
|
|
||
|
|
|
||
|
|
function [model,keptvariables,newmonomtable] = setupQuadratics(model,keptvariables,newmonomtable);
|
||
|
|
if any(model.variabletype) & all(model.variabletype <= 2)
|
||
|
|
monomials = find(model.variabletype);
|
||
|
|
if nnz(model.F_struc(:,1+monomials))==0
|
||
|
|
if all(isinf(model.lb(monomials)))
|
||
|
|
if all(isinf(model.ub(monomials)))
|
||
|
|
% OK, the quadratic/bilinear terms only enter in objective,
|
||
|
|
% so let us try to construct Q
|
||
|
|
if isempty(model.precalc.Qmap)
|
||
|
|
ii = [];
|
||
|
|
jj = [];
|
||
|
|
kk = [];
|
||
|
|
for k = monomials(:)'
|
||
|
|
i = find(model.monomtable(k,:));
|
||
|
|
if model.variabletype(k)==1
|
||
|
|
model.Q(i(1),i(2)) = model.Q(i(1),i(2)) + model.c(k)/2;
|
||
|
|
model.Q(i(2),i(1)) = model.Q(i(1),i(2));
|
||
|
|
ii = [ii;i(1)];
|
||
|
|
jj = [jj;i(2)];
|
||
|
|
kk = [kk;k];
|
||
|
|
else
|
||
|
|
model.Q(i,i) = model.Q(i,i) + model.c(k);
|
||
|
|
ii = [ii;i];
|
||
|
|
jj = [jj;i];
|
||
|
|
kk = [kk;k];
|
||
|
|
end
|
||
|
|
end
|
||
|
|
model.precalc.Qmap.M = sparse(sub2ind(size(model.Q),ii,jj),kk,1,prod(size(model.Q)),length(model.c));
|
||
|
|
model.precalc.Qmap.monomials = monomials;
|
||
|
|
model.precalc.Qmap.monomtable = model.monomtable;
|
||
|
|
else
|
||
|
|
m = model.precalc.Qmap.M*sparse(model.c);
|
||
|
|
m = reshape(m,sqrt(length(m)),[]);
|
||
|
|
model.Q = (m+m')/2;
|
||
|
|
end
|
||
|
|
model.c(monomials)=[];
|
||
|
|
model.F_struc(:,1+monomials) = [];
|
||
|
|
model.lb(monomials) = [];
|
||
|
|
model.ub(monomials) = [];
|
||
|
|
newmonomtable(monomials,:) = [];
|
||
|
|
newmonomtable(:,monomials) = [];
|
||
|
|
model.monomtable = newmonomtable;
|
||
|
|
model.Q(:,monomials) = [];
|
||
|
|
model.Q(monomials,:) = [];
|
||
|
|
model.x0(monomials) = [];
|
||
|
|
model.variabletype(monomials)=[];
|
||
|
|
keptvariables(monomials) = [];
|
||
|
|
end
|
||
|
|
end
|
||
|
|
end
|
||
|
|
end
|