diff --git a/@sdpvar/cos.m b/@sdpvar/cos.m
index 8381321d..a9d99cfa 100644
--- a/@sdpvar/cos.m
+++ b/@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];
diff --git a/@sdpvar/sin.m b/@sdpvar/sin.m
index 5290fe26..54c8a6d6 100644
--- a/@sdpvar/sin.m
+++ b/@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];
diff --git a/@sdpvar/std.m b/@sdpvar/std.m
new file mode 100644
index 00000000..931f573e
--- /dev/null
+++ b/@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
diff --git a/@sdpvar/subsasgn.m b/@sdpvar/subsasgn.m
index 7b4860ce..f5931036 100644
--- a/@sdpvar/subsasgn.m
+++ b/@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
 
diff --git a/@sdpvar/var.m b/@sdpvar/var.m
index 29a7ff47..65ab0f4a 100644
--- a/@sdpvar/var.m
+++ b/@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);
diff --git a/extras/export.m b/extras/export.m
index 3b79a407..de1d92dd 100644
--- a/extras/export.m
+++ b/extras/export.m
@@ -167,6 +167,9 @@
     case 'osqp'
         model = yalmip2osqp(interfacedata);
 
+    case 'piqp'
+        model = yalmip2piqp(interfacedata);
+
     case 'cbc'
         model = yalmip2cbc(interfacedata);
 
diff --git a/extras/sdpsettings.m b/extras/sdpsettings.m
index 7a2ea8dd..5c0aa50b 100644
--- a/extras/sdpsettings.m
+++ b/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;
diff --git a/modules/global/bmibnb.m b/modules/global/bmibnb.m
index 9c19e3f0..b1f13789 100644
--- a/modules/global/bmibnb.m
+++ b/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);
diff --git a/modules/global/compileQuadratic.m b/modules/global/compileQuadratic.m
index 183fbd70..f21e3f51 100644
--- a/modules/global/compileQuadratic.m
+++ b/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
diff --git a/modules/global/preprocess_bilinear_bounds.m b/modules/global/preprocess_bilinear_bounds.m
index e7d3e3c1..007070d4 100644
--- a/modules/global/preprocess_bilinear_bounds.m
+++ b/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 = [];
diff --git a/modules/global/presolve_boundargument_periodic.m b/modules/global/presolve_boundargument_periodic.m
new file mode 100644
index 00000000..ed7eb961
--- /dev/null
+++ b/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
+
+
diff --git a/modules/global/presolve_implied_integer.m b/modules/global/presolve_implied_integer.m
index 6c75c320..36e69d84 100644
--- a/modules/global/presolve_implied_integer.m
+++ b/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);
diff --git a/modules/global/propagate_bounds_from_arbitrary_quadratics.m b/modules/global/propagate_bounds_from_arbitrary_quadratics.m
index d6a1e829..fbb93b8a 100644
--- a/modules/global/propagate_bounds_from_arbitrary_quadratics.m
+++ b/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;
diff --git a/modules/global/update_one_inverseeval_bound.m b/modules/global/update_one_inverseeval_bound.m
index a482000b..df0e4147 100644
--- a/modules/global/update_one_inverseeval_bound.m
+++ b/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
diff --git a/solvers/callpiqp.m b/solvers/callpiqp.m
new file mode 100644
index 00000000..eb7ec75f
--- /dev/null
+++ b/solvers/callpiqp.m
@@ -0,0 +1,69 @@
+function output = callpiqp(interfacedata)
+
+options = interfacedata.options;
+model = yalmip2piqp(interfacedata);
+
+if options.savedebug
+    save debugfile model
+end
+
+if options.showprogress;showprogress(['Calling ' interfacedata.solver.tag],options.showprogress);end
+
+% Define verbose option
+options.piqp.verbose = options.verbose;
+
+% Determine whether problem is sparse
+is_sparse = issparse(model.P) || issparse(model.A) || issparse(model.G);
+
+% Solve with PIQP
+if is_sparse
+    PIQPSolver = piqp('sparse');
+else
+    PIQPSolver = piqp('dense');
+end
+PIQPSolver.setup(model.P, model.c, ...
+                 model.A, model.b, ...
+                 model.G, model.h, ...
+                 model.x_lb, model.x_ub, ...
+                 options.piqp);
+results = PIQPSolver.solve();
+
+switch results.info.status_val
+    case 1
+        problem = 0;
+    case -1
+        problem = 3;
+    case -2
+        problem = 1;
+    case -3
+        problem = 1;
+    case -8
+        problem = 4;
+    case -9
+        problem = 7;
+    case -10
+        problem = 7;
+    otherwise
+        problem = -1;
+end
+
+% Solver time
+solvertime = results.info.run_time;
+
+% Standard interface
+Primal      = results.x(:);
+Dual        = [results.y; results.z];
+infostr     = yalmiperror(problem,interfacedata.solver.tag);
+if ~options.savesolverinput
+    solverinput = [];
+else
+    solverinput = model;
+end
+if ~options.savesolveroutput
+    solveroutput = [];
+else
+    solveroutput = results;
+end
+
+% Standard interface
+output = createOutputStructure(Primal,Dual,[],problem,infostr,solverinput,solveroutput,solvertime);
diff --git a/solvers/definesolvers.m b/solvers/definesolvers.m
index ce6a6703..4332ef88 100644
--- a/solvers/definesolvers.m
+++ b/solvers/definesolvers.m
@@ -710,6 +710,13 @@
 solver(i).constraint.binary = 1;
 i = i+1;
 
+solver(i) = qpsolver;
+solver(i).tag     = 'PIQP';
+solver(i).version = '';
+solver(i).checkfor= {'piqp'};
+solver(i).call    = 'callpiqp';
+i = i+1;
+
 solver(i) = qpsolver;
 solver(i).tag     = 'qpOASES';
 solver(i).version = '';
diff --git a/solvers/yalmip2piqp.m b/solvers/yalmip2piqp.m
new file mode 100644
index 00000000..d9bf1b03
--- /dev/null
+++ b/solvers/yalmip2piqp.m
@@ -0,0 +1,13 @@
+function model = yalmip2piqp(interfacedata);
+
+quadprog_model = yalmip2quadprog(interfacedata);
+
+model.options = interfacedata.options.piqp;
+model.P = quadprog_model.Q;
+model.c = quadprog_model.c;
+model.A = quadprog_model.Aeq;
+model.b = quadprog_model.beq;
+model.G = quadprog_model.A;
+model.h = quadprog_model.b;
+model.x_lb = quadprog_model.lb;
+model.x_ub = quadprog_model.ub;
\ No newline at end of file