Merge branch 'develop' of https://github.com/yalmip/YALMIP into develop

@sdpvar/cos.m

@@ -14,7 +14,7 @@
         operator.monotonicity  = @monotonicity;
         operator.stationary = @stationary; 
         operator.inflection = @inflection;
-        operator.periodic = pi;         
+        operator.periodic = 2*pi;         
         operator.derivative = @(x)(-sin(x));
         operator.range = [-1 1];        
 
@@ -87,6 +87,10 @@
 end
 
 function inflections = inflection(xL,xU)
+if isinf(xL) || isinf(xU)
+    inflections = [];
+    return
+end
 r = floor(xU/(pi/2));
 t = ceil(xL/(pi/2));
 spots = [t:r];
@@ -103,6 +107,11 @@
 end
 
 function [xS,fS] = stationary(xL,xU)
+if isinf(xL) || isinf(xU)
+    xS = [];
+    fS = [];
+    return
+end
 r = floor(xU/(pi/2));
 t = ceil(xL/(pi/2));
 spots = [t:r];

@sdpvar/sin.m

@@ -15,7 +15,7 @@
         operator.bounds     = @bounds;     
         operator.stationary = @stationary; 
         operator.inflection = @inflection;
-        operator.periodic = pi; 
+        operator.periodic = 2*pi; 
         operator.derivative = @(x)(cos(x));
         operator.range = [-1 1];
 
@@ -76,6 +76,10 @@
 end
 
 function inflections = inflection(xL,xU)
+if isinf(xL) || isinf(xU)
+    inflections = [];
+    return
+end
 r = floor(xU/(pi/2));
 t = ceil(xL/(pi/2));
 spots = [t:r];
@@ -91,6 +95,11 @@
     inflections = reshape([xS;dir],1,[]);
 end
 function [xS,fS] = stationary(xL,xU)
+if isinf(xL) || isinf(xU)
+    xS = [];
+    fS = [];
+    return
+end
 r = floor(xU/(pi/2));
 t = ceil(xL/(pi/2));
 spots = [t:r];

@sdpvar/std.m

@@ -0,0 +1,37 @@
+function varargout=std(varargin)
+%STD (overloaded)
+%
+% V = std(x)
+
+switch class(varargin{1})
+
+    case 'sdpvar' % Overloaded operator for SDPVAR objects. Pass on args and save them.
+        X = varargin{1};
+        if nargin > 1 || min(size(X))>1
+            error('SDPVAR/STD only supports simple 1-D variance'),
+        end
+        varargout{1} = yalmip('define',mfilename,varargin{1});
+    case 'char'
+        switch varargin{1}
+            case 'graph'
+                t = varargin{2};
+                X = varargin{3};
+                F = cone([t;X])
+                varargout{1} = F;
+                varargout{2} = struct('convexity','convex','monotonicity','none','definiteness','positive','model','graph');
+                varargout{3} = X;
+
+            case 'exact'
+                t = varargin{2};
+                X = varargin{3};
+                e = sdpvar(length(X),1);
+                F = [e'*e == t^2, e == X, t>=0];
+                varargout{1} = F;
+                varargout{2} = struct('convexity','convex','definiteness','positive','model','exact');
+                varargout{3} = X;
+            otherwise
+                error('SDPVAR/NORM called with weird argument?');
+        end
+    otherwise
+        error('Weird first argument in SDPVAR/STD')
+end

@sdpvar/subsasgn.m

@@ -10,10 +10,14 @@
         end
         if any(I.subs{1} <=0)
             error('Index into matrix is negative or zero.');
-        end 
-
+        end         
         if length(I.subs)>2 && X_is_spdvar
-            y = subsasgn(ndsdpvar(X),I,Y);
+            if isequal(I.subs{end},':')
+                I.subs = {I.subs{1:end-1}};
+                y = subsasgn(X,I,Y);               
+            else
+                y = subsasgn(ndsdpvar(X),I,Y);
+            end           
             return
         end
 

@sdpvar/var.m

@@ -5,13 +5,13 @@
 
 x = varargin{1};
 
-if nargin > 1 | min(size(x))>1
+if nargin > 1 || min(size(x))>1
     error('SDPVAR/VAR only supports simple 1-D variance'),
 end
 
 switch length(x)
     case 1
-        varargout{1} = x;
+        varargout{1} = 0;
     otherwise
         x = reshape(x,length(x),1);
         m = sum(x)/length(x);

extras/export.m

@@ -167,6 +167,9 @@
     case 'osqp'
         model = yalmip2osqp(interfacedata);
 
+    case 'piqp'
+        model = yalmip2piqp(interfacedata);
+
     case 'cbc'
         model = yalmip2cbc(interfacedata);
 

extras/sdpsettings.m

@@ -233,6 +233,9 @@
     options.pensdp = setup_pensdp_options;
     Names = appendOptionNames(Names,options.pensdp,'pensdp');
 
+    options.piqp = setup_piqp_options;
+    Names = appendOptionNames(Names,options.piqp,'piqp');
+
     options.pop = setup_pop_options;
     Names = appendOptionNames(Names,options.pop,'pop');
 
@@ -290,7 +293,7 @@
     options.default.copt = options.copt;
     options.default.cplex = options.cplex;
     options.default.gurobi = options.gurobi;
-    options.default.mosek = options.mosek;   
+    options.default.mosek = options.mosek;
     options.default.osqp = options.osqp;   
     options.default.xpress = options.xpress;   
 end
@@ -1049,6 +1052,14 @@
 pensdp.ALPHA = 1e-2;
 pensdp.P0 = 0.9;
 
+function piqp_options = setup_piqp_options
+try
+    s = piqp('dense');
+    piqp_options = s.get_settings();
+catch
+    piqp_options =[];
+end
+
 function sparsecolo = setup_sparsecolo_options
 sparsecolo.SDPsolver = '';
 sparsecolo.domain = 2;

modules/global/bmibnb.m

@@ -112,6 +112,7 @@
 p.nonshiftedQP.Q =[];
 p.nonshiftedQP.c =[];
 p.nonshiftedQP.f =[];
+p.possibleSol = [];
 % Flags used by BNB which uses common presolve
 p.binarycardinality.up = length(p.binary_variables);
 p.binarycardinality.down = 0;
@@ -203,6 +204,7 @@
 p = detect_hiddendelayedconvex_sdp(p);
 p = presolve_bounds_from_domains(p);
 p = presolve_bounds_from_modelbounds(p,1);
+p = presolve_boundargument_periodic(p);
 
 % *************************************************************************
 % Start some bound propagation
@@ -560,7 +562,7 @@
 else
     if ~isinf(upper)
         if p.options.bmibnb.verbose>0
-            disp('* -Solution proven optimal already in root-node.');   
+            disp('* -Solution proven optimal already during presolve.');   
         end
         counter = p.counter;
         problem = 0;        
@@ -570,6 +572,9 @@
         History.upper = upper;
         History.solutions = x_min;
     else
+        if p.options.bmibnb.verbose>0
+            disp('* -Problem proven infeasible during presolve.');   
+        end
         counter = p.counter;
         problem = 1;
         x_min = repmat(nan,length(p.c),1);

modules/global/compileQuadratic.m

@@ -1,8 +1,8 @@
 function [Q,c] = compileQuadratic(c,p,onlysquares)
 Q = spalloc(length(c),length(c),0);
-if ~p.options.bmibnb.lowerpsdfix
-    return
-end
+% if ~p.options.bmibnb.lowerpsdfix
+%     return
+% end
 for i = 1:size(p.bilinears,1)
     quadratic = (p.bilinears(i,2)==p.bilinears(i,3));
     if onlysquares == 3 && quadratic

modules/global/preprocess_bilinear_bounds.m

@@ -55,10 +55,10 @@
         end
     end
 end
-p.lb(p.binary_variables) = max(0,p.lb(p.binary_variables));
-p.ub(p.binary_variables) = min(1,p.ub(p.binary_variables));
-p.lb(p.integer_variables) = ceil(p.lb(p.integer_variables));
-p.ub(p.integer_variables) = floor(p.ub(p.integer_variables));
+p.lb(p.binary_variables) = max(0,p.lb(p.binary_variables)-1e-6);
+p.ub(p.binary_variables) = min(1,p.ub(p.binary_variables)+1e-6);
+p.lb(p.integer_variables) = ceil(p.lb(p.integer_variables)-1e-6);
+p.ub(p.integer_variables) = floor(p.ub(p.integer_variables)+1e-6);
 p = clean_bounds(p);
 if size(p.F_struc,1)==0
     p.F_struc = [];

modules/global/presolve_boundargument_periodic.m

@@ -0,0 +1,42 @@
+function model = presolve_boundargument_periodic(model)
+
+if ~isempty(model.evalMap)
+    % Find arguments only used in periodic functions, and find largest
+    % period
+    Q = -inf(1,length(model.c));
+    for i = 1:length(model.evalMap)
+        if ~isempty(model.evalMap{i}.properties.periodic)
+            k = model.evalMap{i}.variableIndex;        
+            Q(k) = max(Q(k), model.evalMap{i}.properties.periodic);            
+        else
+            k = model.evalMap{i}.variableIndex;        
+            Q(k) = inf;
+        end
+    end    
+    for i = 1:length(model.evalMap)
+        if ~isempty(model.evalMap{i}.properties.periodic)
+            k = model.evalMap{i}.variableIndex;
+            P = Q(k);
+            if ~isinf(Q(k))                               
+                    if nnz(model.monomtable(:,k)) == 1
+                        % It is  ot used in any monomial
+                        if isempty(model.F_struc) || nnz(model.F_struc(:,k+1))==0
+                            % and not used in any constraint
+                            % Hence, since variable is used nowhere except as
+                            % argument in this periodic function, we can bound
+                            if isinf(model.lb(k)) && isinf(model.ub(k))
+                                model.lb(k) = 0;
+                                model.ub(k) = P;
+                            elseif ~isinf(model.lb(k))
+                                model.ub(k) = min(model.ub(k),model.lb(k)+P);
+                            elseif ~isinf(model.ub(k))
+                                model.lb(k) = max(model.lb(k),model.ub(k)-P);
+                        end
+                    end
+                end
+            end
+        end
+    end
+end
+
+

modules/global/presolve_implied_integer.m

@@ -41,3 +41,5 @@
     end
     p.binary_variables = unique(p.binary_variables );
 end
+p.integer_variables = unique(p.integer_variables );
+p.integer_variables = setdiff(p.integer_variables ,p.binary_variables);

modules/global/propagate_bounds_from_arbitrary_quadratics.m

@@ -35,8 +35,8 @@
                         a(x) = 0;
                         indNEG = find(a < 0);
                         indPOS = find(a > 0);
-                        LB = p.lb;
-                        UB = p.ub;
+                        LB = p.lb-1e-8;
+                        UB = p.ub+1e-8;
                         LB(k) = 0;
                         UB(k) = 0;
                         LB(x) = 0;

modules/global/update_one_inverseeval_bound.m

@@ -9,13 +9,30 @@
     try
         if ~isempty(p.evalMap{i}.properties.inverse)
             if strcmpi(p.evalMap{i}.properties.monotonicity,'increasing')
-                invfiL = real(feval(p.evalMap{i}.properties.inverse,fL,p.evalMap{i}.arg{2:end-1}));
-                invfiU = real(feval(p.evalMap{i}.properties.inverse,fU,p.evalMap{i}.arg{2:end-1}));
+                % Be careful about edge-cases
+                if isequal(fL,p.evalMap{i}.properties.range(1))    
+                    invfiL = p.evalMap{i}.properties.domain(1);
+                else
+                    invfiL = real(feval(p.evalMap{i}.properties.inverse,fL,p.evalMap{i}.arg{2:end-1}));
+                end
+                if isequal(fL,p.evalMap{i}.properties.range(2))    
+                    invfiU = p.evalMap{i}.properties.domain(2);
+                else
+                    invfiU = real(feval(p.evalMap{i}.properties.inverse,fU,p.evalMap{i}.arg{2:end-1}));
+                end
                 p.lb(arg) = max(p.lb(arg),invfiL);
                 p.ub(arg) = min(p.ub(arg),invfiU);
-            elseif strcmpi(p.evalMap{i}.properties.monotonicity,'decreasing')
-                invfiL = real(feval(p.evalMap{i}.properties.inverse,fU,p.evalMap{i}.arg{2:end-1}));
-                invfiU = real(feval(p.evalMap{i}.properties.inverse,fL,p.evalMap{i}.arg{2:end-1}));
+            elseif strcmpi(p.evalMap{i}.properties.monotonicity,'decreasing')                
+                if isequal(fU,p.evalMap{i}.properties.range(2))
+                    invfiL = p.evalMap{i}.properties.domain(1);
+                else
+                    invfiL = real(feval(p.evalMap{i}.properties.inverse,fU,p.evalMap{i}.arg{2:end-1}));
+                end
+                if isequal(fL,p.evalMap{i}.properties.range(1))                    
+                    invfiU = p.evalMap{i}.properties.domain(2);
+                else
+                    invfiU = real(feval(p.evalMap{i}.properties.inverse,fL,p.evalMap{i}.arg{2:end-1}));
+                end
                 p.lb(arg) = max(p.lb(arg),invfiL);
                 p.ub(arg) = min(p.ub(arg),invfiU);
             end