diff --git a/initCobraToolbox.m b/initCobraToolbox.m
index 3404c44ef5..ca323154b3 100644
--- a/initCobraToolbox.m
+++ b/initCobraToolbox.m
@@ -16,11 +16,8 @@ function initCobraToolbox(updateToolbox)
 %
 %     initCobraToolbox
 %           -or-
-%     changeCobraSolver('gurobi');
-%     changeCobraSolver('gurobi', 'MILP');
-%     changeCobraSolver('tomlab_cplex', 'QP');
-%     changeCobraSolver('tomlab_cplex', 'MIQP');
-%     changeCbMapOutput('svg');
+%     changeCobraSolver('gurobi','all');
+%     changeCobraSolver('mosek', 'CLP');
 %
 %     Maintained by Ronan M.T. Fleming
 
@@ -362,18 +359,7 @@ function initCobraToolbox(updateToolbox)
 % requires the solver compatibility to be re-read at each initialisation
 clear isCompatible
 
-%These default tolerances are based on the default values for the Gurobi LP
-%solver. Do not change them without first consulting with other developers.
-%https://www.gurobi.com/documentation/9.0/refman/parameters.html
-% (primal) feasibility tolerance
-changeCobraSolverParams('LP', 'feasTol', 1e-6);
-% (dual) optimality tolerance
-changeCobraSolverParams('LP', 'optTol', 1e-6);
 
-% (primal) feasibility tolerance
-changeCobraSolverParams('EP', 'feasTol', 1e-8);
-% (dual) optimality tolerance
-changeCobraSolverParams('EP', 'optTol', 1e-12);
 
 % Check that SBML toolbox is installed and accessible
 if ~exist('TranslateSBML', 'file')
@@ -409,10 +395,8 @@ function initCobraToolbox(updateToolbox)
 %Define a set of "use first" solvers, other supported solvers will also be added to the struct.
 %This allows to assign them in any order but keep the most commonly used ones on top of the struct.
 SOLVERS = struct('gurobi',struct(),...
-    'ibm_cplex',struct(),...
-    'tomlab_cplex',struct(),...
-    'glpk',struct(),...
     'mosek',struct(),...
+    'glpk',struct(),...
     'matlab',struct());
 
 % active support - supported solvers
@@ -423,17 +407,17 @@ function initCobraToolbox(updateToolbox)
 SOLVERS.quadMinos.type = {'LP'};
 SOLVERS.dqqMinos.type = {'LP','QP'};
 SOLVERS.matlab.type = {'LP', 'NLP'};
-% active support of cplex interfaces - supported solvers
-SOLVERS.cplex_direct.type = {'LP', 'MILP', 'QP'};
-SOLVERS.ibm_cplex.type = {'LP', 'MILP', 'QP', 'MIQP'};
-SOLVERS.cplexlp.type = {'LP'};
-SOLVERS.tomlab_cplex.type = {'LP', 'MILP', 'QP', 'MIQP'};
+
 
 % passive support - solver interfaces
 SOLVERS.qpng.type = {'QP'};
 SOLVERS.tomlab_snopt.type = {'NLP'};
 
 % legacy solvers
+SOLVERS.cplex_direct.type = {'LP', 'MILP', 'QP'};
+SOLVERS.ibm_cplex.type = {'LP', 'MILP', 'QP', 'MIQP'};
+SOLVERS.cplexlp.type = {'LP'};
+SOLVERS.tomlab_cplex.type = {'LP', 'MILP', 'QP', 'MIQP'};
 %SOLVERS.gurobi_mex.type = {'LP', 'MILP', 'QP', 'MIQP'};
 %SOLVERS.lindo_old.type = {'LP'};
 %SOLVERS.lindo_legacy.type = {'LP'};
@@ -441,7 +425,6 @@ function initCobraToolbox(updateToolbox)
 %SOLVERS.opti.type = {'LP', 'MILP', 'QP', 'MIQP', 'NLP'};
 
 % definition of category of solvers with active support
-
 SOLVERS.dqqMinos.categ = 'active';
 SOLVERS.glpk.categ = 'active';
 SOLVERS.gurobi.categ = 'active';
@@ -450,11 +433,6 @@ function initCobraToolbox(updateToolbox)
 SOLVERS.pdco.categ = 'active';
 SOLVERS.quadMinos.categ = 'active';
 
-SOLVERS.cplex_direct.categ = 'active';
-SOLVERS.ibm_cplex.categ = 'active';
-SOLVERS.cplexlp.categ = 'active';
-SOLVERS.tomlab_cplex.categ = 'active';
-
 % definition of category of solvers with passive support
 SOLVERS.qpng.categ = 'passive';
 SOLVERS.tomlab_snopt.categ = 'passive';
@@ -464,6 +442,10 @@ function initCobraToolbox(updateToolbox)
 %SOLVERS.lindo_old.categ = 'legacy';
 %SOLVERS.lindo_legacy.categ = 'legacy';
 SOLVERS.lp_solve.categ = 'legacy';
+SOLVERS.cplex_direct.categ = 'legacy';
+SOLVERS.ibm_cplex.categ = 'legacy';
+SOLVERS.cplexlp.categ = 'legacy';
+SOLVERS.tomlab_cplex.categ = 'legacy';
 %SOLVERS.opti.categ = 'legacy';
 
 % definition of categories of solvers
@@ -484,17 +466,19 @@ function initCobraToolbox(updateToolbox)
     end
 end
 
-% check the installation of the solver
+supportedSolversNames = setdiff(supportedSolversNames,{'cplex_direct','ibm_cplex','cplexlp','tomlab_cplex'});
+
+% check the installation of the solver - except cplex
 for i = 1:length(supportedSolversNames)
     if 0 %set to 1 to debug a new solver
         disp(supportedSolversNames{i})
-        if strcmp(supportedSolversNames{i},'quadMinos')
-            pause(0.1)
+        if strcmp(supportedSolversNames{i},'quadMinos') 
+            disp(supportedSolversNames{i})
         end
     end
     %We will validate all solvers in init. After this, all solvers are
     %checked, whether they actually work and the SOLVERS field is set.
-    [solverOK,solverInstalled] = changeCobraSolver(supportedSolversNames{i},SOLVERS.(supportedSolversNames{i}).type{1},0, 2);
+    [solverOK,solverInstalled] = changeCobraSolver(supportedSolversNames{i},SOLVERS.(supportedSolversNames{i}).type{1},0, 1);
     if strcmp(supportedSolversNames{i},'gurobi') && 0%use fordebugging
         disp(supportedSolversNames{i});
     end
@@ -667,7 +651,6 @@ function initCobraToolbox(updateToolbox)
 
 % use Gurobi (if installed) as the default solver for LP, QP and MILP problems
 changeCobraSolver('gurobi', 'ALL', 0);
-%changeCobraSolver('ibm_cplex', 'QP', 0); %until problem with gurobi QP sorted
 
 % check if a new update exists
 if installedGit && ENV_VARS.printLevel && status_curl == 0 && contains(result_curl, ' 200') && updateToolbox
diff --git a/src/analysis/FBA/optimizeCbModel.m b/src/analysis/FBA/optimizeCbModel.m
index 731bcdb1a9..72a9b8bdcf 100644
--- a/src/analysis/FBA/optimizeCbModel.m
+++ b/src/analysis/FBA/optimizeCbModel.m
@@ -61,6 +61,13 @@
 %                         * ctrs `k x 1` Cell Array of Strings giving IDs of the coupling constraints
 %                         * dsense - `k x 1` character array with entries in {L,E,G}
 %
+%                         * `.evars` : evars x 1  Column Cell Array of Strings	IDs of the additional variables
+%                         * `.E`     : n x evars  The additional Variable Matrix
+%                         * `.evarub`: evars x 1  The upper bounds of the variables from E;
+%                         * `.evarlb`: evars x 1  The lower bounds of the variables from E;
+%                         * `.evarc` : evars x 1  The objective coefficients of the variables from E;
+%                         * `.D`     : k x evars  The matrix coupling additional Constraints (form C), with additional Variables (from E);
+%
 %                         * g0 - `n x 1` weights on zero norm, where positive is minimisation, negative is maximisation, zero is neither.
 %                         * g1 - `n x 1` weights on one norm, where positive is minimisation, negative is maximisation, zero is neither.
 %                         * g2 - `n x 1` weights on two norm
@@ -152,11 +159,11 @@
 %
 % OUTPUT:
 %    solution:       solution object:
-%
+%                          * obj - Primal objective value (sum of f0,f1,f2 terms, ignoring NaN)
 %                          * f - Linear objective value (from LP problem)
 %                          * f0 - Zero-norm objective value
-%                          * f1 - One-norm objective value 
-%                          * f2 - Two-norm objective value
+%                          * f1 - Linear part of objective value (c'*v or ||v||_1)
+%                          * f2 - Quadratic part of objective value
 %                          * v - Reaction rates (Optimal primal variable, legacy FBAsolution.x)
 %                          * y - Dual to the matrix inequality constraints (Shadow prices)
 %                          * w - Dual to the box constraints (Reduced costs)
@@ -222,18 +229,41 @@
 
 
 % Figure out linear objective sense
-if exist('osenseStr', 'var') 
+if exist('osenseStr', 'var')
     if isempty(osenseStr)
         model.osenseStr = 'max';
     else
-        model.osenseStr = osenseStr;
-    end
-else
-    if isfield(model, 'osenseStr')
-        model.osenseStr = model.osenseStr;
-    else
-        model.osenseStr = 'max';
+        %second argument may be the parameter structure
+        if isstruct(osenseStr)
+            if exist('param','var')
+                error('osenseStr is a structure and param structure is also present')
+            else
+                %second argument is be the parameter structure
+                param = osenseStr;
+                if isfield(param,'osenseStr') && ~isfield(model, 'osenseStr')
+                    model.osenseStr = param.osenseStr;
+                else
+                    model.osenseStr = 'max';
+                end
+                if isfield(param,'minNorm')
+                    minNorm=param.minNorm;
+                else
+                    minNorm=[];
+                end
+                if isfield(param,'allowLoops')
+                    allowLoops=param.allowLoops;
+                else
+                    allowLoops=1;
+                end
+            end
+        end
     end
+    % % override if osenseStr already in the model
+    % if isfield(model, 'osenseStr')
+    %     model.osenseStr = model.osenseStr;
+    % else
+    %     model.osenseStr = 'max';
+    % end
 end
 
 if ~exist('param','var')
@@ -246,7 +276,7 @@
         if minNorm == true
             minNorm = 1e-6;
         else
-            minNorm = 0;
+            minNorm = [];
         end
     end
     
@@ -367,31 +397,31 @@
 end
 
 %by default, do linear optimisation unless not required
-doLinearOptimisation = 1;
+doLinearOptimisationFirst = 1;
 %if there is no linear objective, do a QP
 if all(model.c==0)
-    doLinearOptimisation = 0;
+    doLinearOptimisationFirst = 0;
 end
 
 % if there is no linear objective and no quadratic objective, do an LP
 if isempty(minNorm)
-    doLinearOptimisation = 1;
+    doLinearOptimisationFirst = 1;
 end
 
 % If there is are linear and quadrative objectives but the bounds on the
 % corresponding reaction are fixed, then there is no need to solve an LP
 % first, so do a QP
 if all( (model.lb == model.ub & model.c~=0) == (model.c~=0)) && ~isempty(minNorm)
-    doLinearOptimisation = 0;
+    doLinearOptimisationFirst = 0;
 end
 
 % If this is a quadratically regularised LP, go straight to QP
 % TODO This is a hack of the param.minNorm to direct solution to QRLP or QRQP
 if isfield(param,'solveWBMmethod')
-    if any(strcmp(param.solveWBMmethod,{'QRLP','QRQP'}))
-        doLinearOptimisation = 0;
-        param.minNormWBM = minNorm;
-        minNorm = param.solveWBMmethod;
+    if any(strcmp(param.solveWBMmethod,{'QP','QRLP','QRQP'}))%TODO  'QRLP','QRQP' need coded in
+        model.c(:)=0;
+        doLinearOptimisationFirst = 0;
+        minNorm = param.minNorm;
     else
         param.solveWBMmethod = [];
     end
@@ -401,13 +431,15 @@
 
 
 % build the optimization problem
-optProblem = buildOptProblemFromModel(model,verify);
+optProblem = buildOptProblemFromModel(model,verify,param);
 % save the original size of the problem
 [~,nTotalVars] = size(optProblem.A); % nTotalVars needed even if optProblem not used for an LP
 
+
+
 %%
 t1 = clock;
-if doLinearOptimisation
+if doLinearOptimisationFirst
     if allowLoops && ~strcmp(minNorm,'optimizeCardinality')
         clear model
     end
@@ -444,7 +476,7 @@
 
 %only run if minNorm is not empty, and either there is no linear objective
 %or there is a linear objective and the LP problem solved to optimality
-if (doLinearOptimisation==0 && ~isempty(minNorm)) || (doLinearOptimisation==1 && solution.stat==1 && ~isempty(minNorm))
+if (doLinearOptimisationFirst==0 && ~isempty(minNorm)) || (doLinearOptimisationFirst==1 && solution.stat==1 && ~isempty(minNorm))
 
     if strcmp(minNorm, 'optimizeCardinality')
         % DC programming for solving the cardinality optimization problem
@@ -519,7 +551,7 @@
         %                   * .thetaMultiplier - at each iteration: theta = theta*thetaMultiplier
         %                   * .eta - Smallest value considered non-zero (Default value feasTol*1000)
         
-        if doLinearOptimisation
+        if doLinearOptimisationFirst
             optProblem2 = optProblem;
             optProblem2.A = [optProblem.A ; optProblem.c'];
             optProblem2.b = [optProblem.b ; objectiveLP];
@@ -555,7 +587,7 @@
         %                   lb <= v <= ub
         
         % Define the constraints structure
-        if doLinearOptimisation
+        if doLinearOptimisationFirst
             optProblem2.A = [optProblem.A ; optProblem.c'];
             optProblem2.b = [optProblem.b ; objectiveLP];
             optProblem2.csense = [optProblem.csense;'E'];
@@ -676,8 +708,8 @@
         solution = solveCobraLP(optProblem2, param);
 
     elseif strcmp(minNorm, 'QRLP')
-
-        buildOptProblemFromModel_param.method='QRLP';
+        buildOptProblemFromModel_param = param;
+        buildOptProblemFromModel_param.minNorm = param.minNormWBM;
         optProblem = buildOptProblemFromModel(model, 0, buildOptProblemFromModel_param);
         solutionQRLP = solveCobraQP(optProblem,param);
 
@@ -698,9 +730,8 @@
         solution.q = -solution.q; % lb + p <= x <= ub + q
 
     elseif strcmp(minNorm, 'QRQP')
-
-        buildOptProblemFromModel_param.method='QRQP';
-        optProblem = buildOptProblemFromModel(model, 0, buildOptProblemFromModel_param);
+        minNorm = minNormWBM;
+        optProblem = buildOptProblemFromModel(model, 0, param);
         solution = solveCobraQP(optProblem);
 
     elseif length(minNorm)> 1 || minNorm > 0
@@ -737,7 +768,7 @@
         end
 
         % quadratic minimization of the norm.
-        if doLinearOptimisation
+        if doLinearOptimisationFirst
             % set previous linear optimum as constraint.
             optProblem2 = optProblem;
             optProblem2.A = [optProblem.A;optProblem.c'];
@@ -749,7 +780,7 @@
                 %quadratic optimization will get rid of the loops unless you are maximizing a flux which is
                 %part of a loop. By definition, exchange reactions are not part of these loops, more
                 %properly called stoichiometrically balanced cycles.
-                solution = solveCobraQP(optProblem2);
+                solution = solveCobraQP(optProblem2,param);
             else
                 %this is slow, but more useful than minimizing the Euclidean norm if one is trying to
                 %maximize the flux through a reaction in a loop. e.g. in flux variablity analysis
@@ -757,14 +788,13 @@
                 solution = solveCobraMIQP(MIQPproblem);
             end
         else
-
             optProblem.F = spdiags(minNorm,0,nTotalVars,nTotalVars);
             if allowLoops
                 %quadratic optimization will get rid of the loops unless you are maximizing a flux which is
                 %part of a loop. By definition, exchange reactions are not part of these loops, more
                 %properly called stoichiometrically balanced cycles.
 
-                solution = solveCobraQP(optProblem);
+                solution = solveCobraQP(optProblem,param);
             else
                 %this is slow, but more useful than minimizing the Euclidean norm if one is trying to
                 %maximize the flux through a reaction in a loop. e.g. in flux variablity analysis
@@ -775,10 +805,10 @@
     end
 end
 
-%TODO fix this Hack in case param.minNorm is used again
-if ~isempty(param.solveWBMmethod)
-    minNorm = param.minNormWBM;
-end
+% %TODO fix this Hack in case param.minNorm is used again
+% if ~isempty(param.solveWBMmethod)
+%     param.minNorm = param.minNormWBM;
+% end
 
 %dummy parts of the solution
 solution.f0 = NaN;
@@ -835,23 +865,35 @@
     %the value of the second part of the objective depends on the norm
     switch minNorm
         case 'empty'
+            solution.f0 = NaN;
             solution.f1 = optProblem.c'*solution.full(1:nTotalVars,1);
+            solution.f2 = NaN;
         case 'zero'
             %zero norm
             solution.f0 = sum(abs(solution.full(1:nTotalVars,1)) > feasTol);
+            solution.f1 = NaN;
+            solution.f2 = NaN;
         case 'one'
+            solution.f0 = NaN;
             %one norm
             solution.f1 = sum(abs(solution.full(1:nTotalVars,1)));
+            solution.f2 = NaN;
         case 'two'
+            solution.f0 = NaN;
             if isfield(optProblem,'c')
                 solution.f1 = optProblem.c'*solution.full(1:nTotalVars,1);
                 if isfield(solution,'objLinear')
                     solution = rmfield(solution,'objLinear');
                 end
             else
-                solution.f1 = 0;
+                solution.f1 = NaN;
+            end
+            if isfield(optProblem,'F')
+                solution.f2 = 0.5*solution.full'*optProblem.F*solution.full;
+            else
+                disp(param)
+                warning('optProblem.F missing')
             end
-            solution.f2 = 0.5*solution.full'*optProblem.F*solution.full;
             if isfield(solution,'objQuadratic')
                 solution.f2 = solution.objQuadratic;
                 solution = rmfield(solution,'objQuadratic');
@@ -859,12 +901,17 @@
         otherwise
             if exist('LPproblem2','var')
                 if isfield(optProblem2,'F')
-                    solution.f0 = 0;
+                    solution.f0 = NaN;
                     solution.f1 = optProblem.c'*solution.full(1:nTotalVars,1);
                     solution.f2 = 0.5*solution.full'*optProblem2.F*solution.full;
                 end
+            else
+                solution.f0 = NaN;
+                solution.f1 = NaN;
+                solution.f2 = NaN;
             end
     end
+    solution.obj = sum([solution.f0,solution.f1,solution.f2],'omitnan');
     
     %primal optimal variables
     solution.v = solution.full(1:nRxns);
@@ -922,7 +969,7 @@
     
     solution.time = etime(clock, t1);
     
-    fieldOrder = {'f';'f0';'f1';'f2';'v';'y';'w';'s';'solver';'method';'stat';'origStat';'time';'basis';'vars_v';'vars_w';'ctrs_y';'ctrs_s';'x';'full';'obj';'rcost';'dual';'slack'};
+    fieldOrder = {'f';'f0';'f1';'f2';'v';'y';'w';'s';'solver';'lpmethod';'qpmethod';'stat';'origStat';'time';'basis';'vars_v';'vars_w';'ctrs_y';'ctrs_s';'x';'full';'obj';'rcost';'dual';'slack'};
     % reorder fields for better readability
     currentfields = fieldnames(solution);
     presentfields = ismember(fieldOrder,currentfields);
@@ -932,6 +979,7 @@
     if 0
         %return NaN of correct dimensions if problem does not solve properly
         solution.f = NaN;
+
         solution.v = NaN*ones(nRxns,1);
         solution.y = NaN*ones(nMets,1);
         solution.w = NaN*ones(nRxns,1);
@@ -947,7 +995,11 @@
     else
         %return empty fields if problem does not solve properly (backward
         %compatible)
+        solution.obj = NaN;
         solution.f = NaN;
+        solution.f0 = NaN;
+        solution.f1 = NaN;
+        solution.f2 = NaN;
         solution.v = [];
         solution.y = [];
         solution.w = [];
@@ -968,9 +1020,6 @@
 if 1 %this may not be very backward compatible
     %remove fields coming from solveCobraLP/QP but not part of the specification
     %of the output from optimizeCbModel
-    if isfield(solution,'obj')
-        solution = rmfield(solution,'obj');
-    end
     if isfield(solution,'full')
         solution = rmfield(solution,'full');
     end
diff --git a/src/analysis/wholeBody/PSCMToolbox/optimizeWBModel.m b/src/analysis/wholeBody/PSCMToolbox/optimizeWBModel.m
index e9e788ae7d..107f4445c7 100644
--- a/src/analysis/wholeBody/PSCMToolbox/optimizeWBModel.m
+++ b/src/analysis/wholeBody/PSCMToolbox/optimizeWBModel.m
@@ -177,7 +177,7 @@
         % is especially well suited for models with a wide range of coefficients in the constraint matrix rows or columns.
         % Settings 1 and 3 are not as directly connected to any specific model characteristics, so experimentation with both
         % settings may be needed to assess performance impact.
-        param.scaleFlag=0;
+        param.ScaleFlag=0;
 
     case 'ibm_cplex'
         % https://www.ibm.com/docs/en/icos/12.10.0?topic=infeasibility-coping-ill-conditioned-problem-handling-unscaled-infeasibilities
diff --git a/src/base/install/updateCobraToolbox.m b/src/base/install/updateCobraToolbox.m
index 5af2667281..cce06c1cb4 100644
--- a/src/base/install/updateCobraToolbox.m
+++ b/src/base/install/updateCobraToolbox.m
@@ -179,7 +179,7 @@ function updateCobraToolbox(fetchAndCheckOnly)
             %disp([codeBaseDir filesep forkNames{j}])
             %disp([toolboxDir filesep dirNames{j}])
             fprintf('%s\n',['> Adding path to ' codeBaseDir filesep forkNames{j} ', and removing path to ' toolboxDir filesep dirNames{j}]);
-            rmpath([toolboxDir filesep dirNames{j}])
+            rmpath(genpath([toolboxDir filesep dirNames{j}]))
             addpath(genpath([codeBaseDir filesep forkNames{j}]))
         end
     end
diff --git a/src/base/solvers/buildOptProblemFromModel.m b/src/base/solvers/buildOptProblemFromModel.m
index 9b88a68365..4ec8f19b31 100644
--- a/src/base/solvers/buildOptProblemFromModel.m
+++ b/src/base/solvers/buildOptProblemFromModel.m
@@ -72,6 +72,8 @@
 %                * `.osense`: Objective sense (`-1`: maximise (default); `1`: minimise)
 %                * `.csense`: string with the constraint sense for each row in A ('E', equality, 'G' greater than, 'L' less than).
 %                * `.F`: Positive semidefinite matrix for quadratic part of objective
+% OPTIONAL OUTPUT:
+
 
 [nMet,nRxn]=size(model.S);
 
@@ -155,13 +157,6 @@
     model.c = zeros(nRxn,1);
 end
 
-if ~isfield(param, 'method')
-    if isfield(model,'F')
-        param.method = 'QP';
-    else
-        param.method = 'LP';
-    end
-end
 
 % case 'LP'
 if ~modelC && ~modelE
@@ -195,86 +190,128 @@
     end
 end
 
-switch param.method
-    case 'LP'
-        %nothing to do - done above already
-    case 'QP'
-        if modelE
-            optProblem.F = spdiags(zeros(size(optProblem.A,2),1),0,size(optProblem.A,2),size(optProblem.A,2));
-            %assume that the remainder of the variables are not being quadratically
-            %minimised
-            optProblem.F(1:size(model.F,1),1:size(model.F,1)) = model.F;
-        else
-            optProblem.F = model.F;
-        end
-    case 'QRLP'
-        [m,n]=size(optProblem.A);
-        optProblem.A = [...
-            %          v           r             p            z
-            optProblem.A,  speye(m,m),  sparse(m,n), sparse(m,n);
-              speye(n,n), sparse(n,m),  -speye(n,n), -speye(n,n)]; % x - p = z;
+if isfield(model,'F')
+    if modelE
+        optProblem.F = spdiags(zeros(size(optProblem.A,2),1),0,size(optProblem.A,2),size(optProblem.A,2));
+        %assume that the remainder of the variables are not being quadratically
+        %minimised
+        optProblem.F(1:size(model.F,1),1:size(model.F,1)) = model.F;
+    else
+        optProblem.F = model.F;
+    end
+end
 
-        optProblem.b  = [optProblem.b;sparse(n,1)];
-        optProblem.lb = [-inf(2*n+m,1); model.lb];
-        optProblem.ub = [ inf(2*n+m,1); model.ub];
-        optProblem.csense = [optProblem.csense;repmat('E',n,1)];
-        optProblem.c = [model.c; sparse(m+2*n,1)];
-        optProblem.F = speye(size(optProblem.A,2),size(optProblem.A,2))*1e-16;%small amound of regularistion makes matrix positive definite numerically
-        optProblem.F(n+1:2*n+m,n+1:2*n+m) = speye(n+m,n+m)*max(abs(model.c))*100; %regularisation must dominate linear objective
-        %dimensions needed to extract non-regularised part of solution
-        optProblem.m = m;
-        optProblem.n = n;
+if isfield(param,'solveWBMmethod') && ~isempty(param.solveWBMmethod)
+    switch param.solveWBMmethod
+        case 'LP'
+            %nothing to do - done above already
+        case 'QP'
+            %nothing to do - done above already
+        case 'QRLP'
+            [m,n]=size(optProblem.A);
+            optProblem.A = [...
+                %          v           r             p            z
+                optProblem.A,  speye(m,m),  sparse(m,n), sparse(m,n);
+                speye(n,n), sparse(n,m),  -speye(n,n), -speye(n,n)]; % x - p = z;
 
-        if 0 %use to test if regularisation enough
-            F2 = optProblem.F;
-            %This line modifies the diagonal elements of F2 by setting them to zero.
-            F2(1:size(optProblem.F,1):end)=0;
-            if all(all(F2)) == 0
-                %only nonzeros in QPproblem.F are on the diagonal
-                try
-                    %try cholesky decomposition
-                    B = chol(optProblem.F);
-                catch
-                    optProblem.F = optProblem.F + diag((diag(optProblem.F)==0)*1e-16);
-                end
-                try
-                    B = chol(optProblem.F);
-                catch
-                    error('QPproblem.F only has non-zeros along the main diagnoal and is still not positive semidefinite after adding 1e-16')
+            optProblem.b  = [optProblem.b;sparse(n,1)];
+            optProblem.lb = [-inf(2*n+m,1); model.lb];
+            optProblem.ub = [ inf(2*n+m,1); model.ub];
+            optProblem.csense = [optProblem.csense;repmat('E',n,1)];
+            optProblem.c = [model.c; sparse(m+2*n,1)];
+            optProblem.F = speye(size(optProblem.A,2),size(optProblem.A,2))*1e-16;%small amound of regularistion makes matrix positive definite numerically
+            optProblem.F(n+1:2*n+m,n+1:2*n+m) = speye(n+m,n+m)*max(abs(model.c))*100; %regularisation must dominate linear objective
+            %dimensions needed to extract non-regularised part of solution
+            optProblem.m = m;
+            optProblem.n = n;
+
+            if 0 %use to test if regularisation enough
+                F2 = optProblem.F;
+                %This line modifies the diagonal elements of F2 by setting them to zero.
+                F2(1:size(optProblem.F,1):end)=0;
+                if all(all(F2)) == 0
+                    %only nonzeros in QPproblem.F are on the diagonal
+                    try
+                        %try cholesky decomposition
+                        B = chol(optProblem.F);
+                    catch
+                        optProblem.F = optProblem.F + diag((diag(optProblem.F)==0)*1e-16);
+                    end
+                    try
+                        B = chol(optProblem.F);
+                    catch
+                        error('QPproblem.F only has non-zeros along the main diagnoal and is still not positive semidefinite after adding 1e-16')
+                    end
                 end
             end
-        end
 
-    case 'QRQP'
-        if modelE
-            optProblem.F = spdiags(zeros(size(optProblem.A,2),1),0,size(optProblem.A,2),size(optProblem.A,2));
-            %assume that the remainder of the variables are not being quadratically
-            %minimised
-            optProblem.F(1:size(model.F,1),1:size(model.F,1)) = model.F;
-        else
-            optProblem.F = model.F;
-        end
+        case 'QRQP'
+            if modelE
+                optProblem.F = spdiags(zeros(size(optProblem.A,2),1),0,size(optProblem.A,2),size(optProblem.A,2));
+                %assume that the remainder of the variables are not being quadratically
+                %minimised
+                optProblem.F(1:size(model.F,1),1:size(model.F,1)) = model.F;
+            else
+                optProblem.F = model.F;
+            end
 
-        [m,n]=size(optProblem.A);
-        optProblem.A = [...
-            optProblem.A,  speye(m,m),  sparse(m,n), sparse(m,n); % A*x + r <=> b
-              speye(n,n), sparse(n,m),  -speye(n,n), -speye(n,n)]; % x - p = z;
+            [m,n]=size(optProblem.A);
+            optProblem.A = [...
+                optProblem.A,  speye(m,m),  sparse(m,n), sparse(m,n); % A*x + r <=> b
+                speye(n,n), sparse(n,m),  -speye(n,n), -speye(n,n)]; % x - p = z;
 
-        optProblem.b  = [optProblem.b;sparse(n,1)];
-        optProblem.lb = [-inf(2*n+m,1); model.lb]; % lb <= x - p       ----> lb <= z
-        optProblem.ub = [ inf(2*n+m,1); model.ub]; %       x - p <= ub ---->       z <= ub
-        optProblem.csense = model.csense;
-        optProblem.c = [model.c; sparse(m+2*n,1)];
-        optProblem.F = sparse(size(optProblem.A,2),size(optProblem.A,2));
-        optProblem.F(1:n,1:n) = model.F;
-        optProblem.F(n+1:2*n+m,n+1:2*n+m) = speye(n+m,n+m)*max(abs(model.c))*100;  %regularisation must dominate linear and quadratic objective
-        %dimensions needed to extract non-regularised part of solution
-        optProblem.m = m;
-        optProblem.n = n;
+            optProblem.b  = [optProblem.b;sparse(n,1)];
+            optProblem.lb = [-inf(2*n+m,1); model.lb]; % lb <= x - p       ----> lb <= z
+            optProblem.ub = [ inf(2*n+m,1); model.ub]; %       x - p <= ub ---->       z <= ub
+            optProblem.csense = model.csense;
+            optProblem.c = [model.c; sparse(m+2*n,1)];
+            optProblem.F = sparse(size(optProblem.A,2),size(optProblem.A,2));
+            optProblem.F(1:n,1:n) = model.F;
+            optProblem.F(n+1:2*n+m,n+1:2*n+m) = speye(n+m,n+m)*max(abs(model.c))*100;  %regularisation must dominate linear and quadratic objective
+            %dimensions needed to extract non-regularised part of solution
+            optProblem.m = m;
+            optProblem.n = n;
 
-    otherwise
-        error('param.method unrecognised')
+        otherwise
+            error('param.method unrecognised')
+    end
 end
 
-
 [~,optProblem.osense] = getObjectiveSense(model);
+
+if isfield(param,'debug') && param.debug
+    switch param.solver
+        case 'mosek'
+            % names
+            % This structure is used to store all the names of individual items in the optimization problem such as the constraints and the variables.
+            %
+            % Fields
+            % name (string) – contains the problem name.
+            %
+            % obj (string) – contains the name of the objective.
+            %
+            % con (cell) – a cell array where names.con{i} contains the name of the
+            % -th constraint.
+            %
+            % var (cell) – a cell array where names.var{j} contains the name of the
+            % -th variable.
+            optProblem.names.name='optimizeCbModel';
+            if any(model.c~=0)
+                optProblem.names.obj=model.rxns{model.c~=0};
+            else
+                optProblem.names.obj='noLP';
+            end
+            if isfield(model,'ctrs')
+                optProblem.names.con=[model.mets;model.ctrs];
+            else
+                optProblem.names.con=model.mets;
+            end
+            if isfield(model,'evars')
+                optProblem.names.var=[model.rxns;model.evars];
+            else
+                optProblem.names.var=model.rxns;
+            end
+    end
+else
+    optProblem.names=[];
+end
\ No newline at end of file
diff --git a/src/base/solvers/cardOpt/sparseLP/optimizeCardinality.m b/src/base/solvers/cardOpt/sparseLP/optimizeCardinality.m
index 45dad396e5..b31014d0bf 100644
--- a/src/base/solvers/cardOpt/sparseLP/optimizeCardinality.m
+++ b/src/base/solvers/cardOpt/sparseLP/optimizeCardinality.m
@@ -91,9 +91,6 @@
 [feasTol] = getCobraSolverParams('LP', 'feasTol');
 
 global CBT_LP_SOLVER
-if strcmp('mosek',CBT_LP_SOLVER)
-    error('mosek not working with optCard  -TODO debug with testOptimizeCbModel')
-end
 
 %% Check inputs
 if ~exist('param','var') || isempty(param)
diff --git a/src/base/solvers/cplex/buildCplexProblemFromCOBRAStruct.m b/src/base/solvers/cplex/buildCplexProblemFromCOBRAStruct.m
index 8891b9bac5..61f1bd766e 100644
--- a/src/base/solvers/cplex/buildCplexProblemFromCOBRAStruct.m
+++ b/src/base/solvers/cplex/buildCplexProblemFromCOBRAStruct.m
@@ -75,7 +75,7 @@
 %always set the problem to minimise, but change the linear objective sign   
 cplexProblem.Model.sense = 'minimize';
 if isfield(Problem,'c')
-    if isfield(Problem,'osense') && isfield(Problem,'c')
+    if isfield(Problem,'osense')
         cplexProblem.Model.obj = Problem.osense*Problem.c;
     else
         cplexProblem.Model.obj = Problem.c;
diff --git a/src/base/solvers/cplex/setCplexParametersForProblem.m b/src/base/solvers/cplex/setCplexParametersForProblem.m
index e4398ab27d..5b1d35b2fa 100644
--- a/src/base/solvers/cplex/setCplexParametersForProblem.m
+++ b/src/base/solvers/cplex/setCplexParametersForProblem.m
@@ -19,6 +19,12 @@
 %set the default parameters so we can see what they are
 cplexProblem.setDefault;
 
+%TODO ? add these 
+% valDef.DATACHECK = 1;
+% valDef.DEPIND = 1;
+% valDef.checkNaN = 0;
+% valDef.warning = 0;
+
 % set the printLevel to the cobra Parameters
 cplexProblem.Param.output.writelevel.Cur = problemTypeParams.printLevel;
 cplexProblem.Param.barrier.display.Cur = problemTypeParams.printLevel;
@@ -67,12 +73,56 @@
 end
 
 % set tolerances
+
+% simplex.tolerances.feasibility
+% Specifies the feasibility tolerance, that is, the degree to which values of the basic variables 
+% calculated by the simplex method may violate their bounds. CPLEX® does not use this tolerance to
+% relax the variable bounds nor to relax right hand side values. This parameter specifies an 
+% allowable violation. Feasibility influences the selection of an optimal basis and can be reset 
+% to a higher value when a problem is having difficulty maintaining feasibility during optimization. 
+% You can also lower this tolerance after finding an optimal solution if there is any doubt 
+% that the solution is truly optimal. If the feasibility tolerance is set too low, CPLEX may falsely
+% conclude that a problem is infeasible. If you encounter reports of infeasibility during Phase II of
+% the optimization, a small adjustment in the feasibility tolerance may improve performance.
+% Values
+% Any number from 1e-9 to 1e-1; default: 1e-06.
 cplexProblem.Param.simplex.tolerances.feasibility.Cur = problemTypeParams.feasTol;
-cplexProblem.Param.simplex.tolerances.optimality.Cur = problemTypeParams.optTol;
 
+
+% network.tolerances.feasibility
+% Specifies feasibility tolerance for network primal optimization. The feasibility tolerance specifies
+% the degree to which the flow value of a model may violate its bounds. This tolerance influences
+% the selection of an optimal basis and can be reset to a higher value when a problem is having 
+% difficulty maintaining feasibility during optimization. You may also wish to lower this tolerance
+% after finding an optimal solution if there is any doubt that the solution is truly optimal. 
+% If the feasibility tolerance is set too low, CPLEX may falsely conclude that a problem is infeasible.
+% If you encounter reports of infeasibility during Phase II of the optimization, a small adjustment
+% in the feasibility tolerance may improve performance.
+% Values
+% Any number from 1e-11 to 1e-1; default: 1e-6.
 cplexProblem.Param.network.tolerances.feasibility.Cur = problemTypeParams.feasTol;
+
+% Influences the reduced-cost tolerance for optimality. This parameter governs 
+% how closely CPLEX must approach the theoretically optimal solution.
+% The simplex algorithm halts when it has found a basic feasible solution with
+% all reduced costs nonnegative. CPLEX uses this optimality tolerance to make 
+% the decision of whether or not a given reduced cost should be considered nonnegative.
+% CPLEX considers "nonnegative" a negative reduced cost having absolute value less
+% than the optimality tolerance. For example, if your optimality tolerance is set
+% to 1e-6, then CPLEX considers a reduced cost of -1e-9 as nonnegative for 
+% the purpose of deciding whether the solution is optimal.
+% Values
+% Any number from 1e-9 to 1e-1; default: 1e-06.
+cplexProblem.Param.simplex.tolerances.optimality.Cur = problemTypeParams.optTol;
+
+% network.tolerances.optimality
+% Specifies the optimality tolerance for network optimization; that is, 
+% the amount a reduced cost may violate the criterion for an optimal solution.
+% Values
+% Any number from 1e-11 to 1e-1; default: 1e-6.
 cplexProblem.Param.network.tolerances.optimality.Cur = problemTypeParams.optTol;
 
+
 %https://www.ibm.com/support/knowledgecenter/SSSA5P_12.7.0/ilog.odms.cplex.help/CPLEX/Parameters/topics/BarEpComp.html
 %Sets the tolerance on complementarity for convergence. The barrier algorithm terminates with an optimal solution if the relative complementarity is smaller than this value.
 %Changing this tolerance to a smaller value may result in greater numerical precision of the solution, but also increases the chance of failure to converge in the algorithm and consequently may result in no solution at all. Therefore, caution is advised in deviating from the default setting.
@@ -215,6 +265,25 @@
     cplexProblem.Param.lpmethod.Cur=4;%BARRIER provided best benchmark performance on Harvetta
 end
 
+if isfield(solverParams,'multiscale') && solverParams.multiscale==1 && 0
+    % Decides how to scale the problem matrix.
+    % Value  Meaning
+    % -1	No scaling
+    % 0	Equilibration scaling; default
+    % 1	More aggressive scaling
+    % https://www.ibm.com/docs/en/icos/12.10.0?topic=parameters-scale-parameter
+    cplexProblem.Param.read.scale.Cur = -1;
+
+    % Emphasizes precision in numerically unstable or difficult problems.
+    % This parameter lets you specify to CPLEX that it should emphasize precision in
+    % numerically difficult or unstable problems, with consequent performance trade-offs in time and memory.
+    % Value Meaning
+    % 0   Do not emphasize numerical precision; default
+    % 1	Exercise extreme caution in computation
+    % https://www.ibm.com/docs/en/icos/12.10.0?topic=parameters-numerical-precision-emphasis
+    cplexProblem.Param.emphasis.numerical.Cur = 1;
+end
+
 if isfield(solverParams,'scaind')
     % Decides how to scale the problem matrix.
     % Value  Meaning
diff --git a/src/base/solvers/entropicFBA/entropicFluxBalanceAnalysis.m b/src/base/solvers/entropicFBA/entropicFluxBalanceAnalysis.m
index 4e74b983b6..8894c2e8bf 100644
--- a/src/base/solvers/entropicFBA/entropicFluxBalanceAnalysis.m
+++ b/src/base/solvers/entropicFBA/entropicFluxBalanceAnalysis.m
@@ -1,5 +1,8 @@
 function [solution, modelOut] = entropicFluxBalanceAnalysis(model, param)
-%% TBC
+% Entropy maximisation of fluxes (or fluxes and concentrations) subject to
+% mass balance, optionally coupling constraints, optionally quadratic
+% penalisation of deviation from given fluxes.
+% 
 % minimize             g.*vf'*(log(vf) -1) + (cf + ci)'*vf 
 % vf,vr,w,x,x0       + g.*vr'*(log(vr) -1) + (cr - ci)'*vr
 %                    + f.*x' *(log(x)  -1) + u0'*x 
@@ -219,6 +222,36 @@
 %% processing for fluxes
 [vl,vu,vel,veu,vfl,vfu,vrl,vru,ci,ce,cf,cr,g] = processFluxConstraints(model,param);
 
+if param.debug && 0 %TODO - remove this
+    modelProcessed = model;
+    modelProcessed.lb(model.SConsistentRxnBool)=vl;
+    modelProcessed.lb(~model.SConsistentRxnBool)=vel;
+    modelProcessed.ub(model.SConsistentRxnBool)=vu;
+    modelProcessed.ub(~model.SConsistentRxnBool)=veu;
+    solutionLP = optimizeCbModel(modelProcessed,param);
+
+    switch solutionLP.stat
+        case 0
+            solution = solutionLP;
+            message = ['solveCobraEP: LP part of EPproblem is infeasible according to solveCobraLP with ' param.solver '.'];
+            warning(message)
+
+            return
+        case 2
+            solution = solutionLP;
+            message = ['solveCobraEP: LP part of EPproblem is unbounded according to solveCobraLP with ' param.solver '.'];
+            warning(message)
+
+            return
+        case 1
+            message =['solveCobraEP: LP part of EPproblem is feasible according to solveCobraLP with ' param.solver '.'];
+            fprintf('%s\n',message)
+        otherwise
+            error('inconclusive solveCobraLP')
+    end
+    messages = cellstr(message);
+
+end
 %% optionally processing for concentrations
 %processConcConstraints
 if contains(lower(param.entropicFBAMethod),'conc')
@@ -1096,25 +1129,11 @@
                 expConeBool = EPproblem.d~=0;
                 nExpCone  = nnz(expConeBool);
                 
-                %
-                if 1
-                    mosekParam=param;
-                    mosekParam.printLevel=param.printLevel-1;
-                    solution = solveCobraEP(EPproblem,mosekParam);
-                else
-                    [verify,method,printLevel,debug,feasTol,optTol,solver,param] =...
-                        getCobraSolverParams('EP',getCobraSolverParamsOptionsForType('EP'),param);
-                    
-                    solution = solveCobraEP(EPproblem,...
-                        'verify',verify,...
-                        'method',method,...
-                        'printLevel',printLevel,...
-                        'debug',debug,...
-                        'feasTol',feasTol,...
-                        'optTol',optTol,...
-                        'solver',solver,...
-                        param);
-                end
+
+                solveCobraEPparam=param;
+                solveCobraEPparam.printLevel=solveCobraEPparam.printLevel-1;
+                solution = solveCobraEP(EPproblem,solveCobraEPparam);
+
                 
                 switch solution.stat
                     case 1
@@ -1606,7 +1625,9 @@
             solution.messages = [];
         end
     otherwise
-        solution_optimizeCbModel = optimizeCbModel(model);
+        %solution = optimizeCbModel(model, osenseStr, minNorm, allowLoops, param)
+        param.debug=1;
+        solution_optimizeCbModel = optimizeCbModel(model,'min',[],1,param);
         switch solution_optimizeCbModel.stat
             case 0
                 message = 'entropicFluxBalanceAnalysis: EPproblem is not feasible, because LP part of model is not feasible according to optimizeCbModel.';
diff --git a/src/base/solvers/entropicFBA/mosekParamSetEFBA.m b/src/base/solvers/entropicFBA/mosekParamSetEFBA.m
deleted file mode 100644
index c435a5708e..0000000000
--- a/src/base/solvers/entropicFBA/mosekParamSetEFBA.m
+++ /dev/null
@@ -1,125 +0,0 @@
-function param=mosekParamSetEFBA(param)
-%creates a structure of pertinent user defined options for MOSEK
-%OUTPUT
-%param      parameter structure to be passed to the MOSEK solver
-
-%MSK_IPAR_LOG_PRESOLVE
-% Description:Controls amount of output printed by the presolve procedure. A higher level implies that more information is logged.
-% Possible Values:Any number between 0 and +inf. 
-% Default value:1
-param.MSK_IPAR_LOG_PRESOLVE =1;
-
-%MSK_IPAR_INTPNT_SCALING 
-% Controls how the problem is scaled before the interior-point optimizer is used.
-% Possible Values:
-%     MSK_SCALING_NONE
-%         No scaling is performed. 
-%     MSK_SCALING_MODERATE
-%         A conservative scaling is performed. 
-%     MSK_SCALING_AGGRESSIVE
-%         A very aggressive scaling is performed. 
-%     MSK_SCALING_FREE
-%         The optimizer chooses the scaling heuristic. 
-% Default value:
-%     MSK_SCALING_FREE
-param.MSK_IPAR_INTPNT_SCALING ='MSK_SCALING_FREE';
-%param.MSK_IPAR_INTPNT_SCALING ='MSK_SCALING_NONE';
-
-% MSK_IPAR_INTPNT_REGULARIZATION_USE 
-% Description:Controls whether regularization is allowed.
-% Possible Values: MSK_ON    Switch the option on. 
-%                  MSK_OFF   Switch the option off. 
-% Default value:   MSK_ON
-param.MSK_IPAR_INTPNT_REGULARIZATION_USE='MSK_ON'; 
-
-%%%%%%%%%%%%% NONLINEAR TERMINATION CRITERIA%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-% MSK_DPAR_INTPNT_CO_TOL_DFEAS 
-%Dual feasibility tolerance used by the interior-point optimizer for conic problems.
-%Default:1.0e-8
-%Accepted: [0.0; 1.0]
-param.MSK_DPAR_INTPNT_CO_TOL_DFEAS = 1e-10;
-
-%MSK_DPAR_INTPNT_CO_TOL_PFEAS
-%Primal feasibility tolerance used by the interior-point optimizer for conic problems.
-%Default: 1.0e-8
-% Accepted: [0.0; 1.0]
-param.MSK_DPAR_INTPNT_CO_TOL_PFEAS = 1.0e-10; %was 1e-11 may be too aggressive -RF
-
-%MSK_DPAR_INTPNT_CO_TOL_REL_GAP
-%Relative gap termination tolerance used by the interior-point optimizer for conic problems.
-%Default:1.0e-8
-%Accepted: [0.0; 1.0]
-param.MSK_DPAR_INTPNT_CO_TOL_REL_GAP = 1.0e-9; %was 1e-11 may be too aggressive -RF
-
-%useful for ensuring dual feasibility is as good as primal
-
-% MSK_IPAR_INTPNT_MAX_ITERATIONS 
-% Controls the maximum number of iterations allowed in the interior-point optimizer.
-% Possible Values:Any number between 0 and +inf. 
-% Default value: 400
-param.MSK_IPAR_INTPNT_MAX_ITERATIONS=400;
-
-
-%%%%%%%%%%%%% NONLINEAR SOLVER INTEGER PARAM %%%%%%%%%%%%%%%%%%%%%%%%%%%%
-
-% MSK_IPAR_BI_IGNORE_MAX_ITER 
-% If the parameter MSK_IPAR_INTPNT_BASIS has the value MSK_BI_NO_ERROR and 
-% the interior-point optimizer has terminated due to maximum number of 
-% iterations, then basis identification is performed if this parameter has
-% the value MSK_ON.
-% Possible Values:
-%     MSK_ON        Switch the option on. 
-%     MSK_OFF       Switch the option off. 
-% Default value:
-%     MSK_OFF 
-param.MSK_IPAR_BI_IGNORE_MAX_ITER='MSK_OFF';     
-
-
-%%%%%%%%%%% Solution Approach
-% MSK_IPAR_INTPNT_SOLVE_FORM 
-% Controls whether the primal or the dual problem is solved.
-% Possible Values:
-%     MSK_SOLVE_PRIMAL
-%         The optimizer should solve the primal problem. 
-%     MSK_SOLVE_DUAL
-%         The optimizer should solve the dual problem. 
-%     MSK_SOLVE_FREE
-%         The optimizer is free to solve either the primal or the dual problem. 
-% Default value:MSK_SOLVE_FREE
-param.MSK_IPAR_INTPNT_SOLVE_FORM='MSK_SOLVE_FREE';
-%param.MSK_IPAR_INTPNT_SOLVE_FORM='MSK_SOLVE_PRIMAL';
-
-%%%%%%% Infeasibility
-% MSK_DPAR_INTPNT_TOL_INFEAS 
-% Controls when the optimizer declares the model primal or dual infeasible.
-% A small number means the optimizer gets more conservative about declaring the model infeasible.
-% Possible Values:Any number between 0.0 and 1.0. 
-% Default value: 1.0e-8
-param.MSK_DPAR_INTPNT_TOL_INFEAS=1e-10; 
-
-%%%%%%%%%%%%%%%%%%%%%%OUTPUT%%%%%%%%%%%%%%%%%%
-%MSK_IPAR_LOG_INTPNT
-% Controls amount of output printed printed by the interior-point optimizer.
-%A higher level implies that more information is logged.
-% Possible Values: Any number between 0 and +inf. 
-% Default value: 4
-param.MSK_IPAR_LOG_INTPNT=5;
-
-%infesibility report
-% MSK_IPAR_INFEAS_REPORT_AUTO
-%Controls the amount of information presented in an infeasibility report.
-% Possible Values:
-%     MSK_ON
-%         Switch the option on. 
-%     MSK_OFF
-%         Switch the option off. 
-% Default value:
-%     MSK_OFF
-param.MSK_IPAR_INFEAS_REPORT_AUTO='MSK_OFF'; 
-
-% MSK_IPAR_INFEAS_REPORT_LEVEL 
-% Controls the amount of information presented in an infeasibility report. Higher values imply more information.
-% Possible Values:Any number between 0 and +inf. 
-% Default value: 1
-%Higher values imply more information.
-param.MSK_IPAR_INFEAS_REPORT_LEVEL=100; 
diff --git a/src/base/solvers/entropicFBA/processConcConstraints.m b/src/base/solvers/entropicFBA/processConcConstraints.m
index 57caf4398c..ae486af86c 100644
--- a/src/base/solvers/entropicFBA/processConcConstraints.m
+++ b/src/base/solvers/entropicFBA/processConcConstraints.m
@@ -27,8 +27,6 @@
 %  model.dcu:     m x 1    real valued upper bound on difference between final and initial initial molecular concentrations  (default inf)
 %  model.gasConstant:    scalar gas constant (default 8.31446261815324 J K^-1 mol^-1)
 %  model.temperature:              scalar temperature (default 310.15 Kelvin)
-%
-%  param.method:  'fluxConc'
 %  param.maxConc: (1e4) maximim micromolar concentration allowed
 %  param.maxConc: (1e-4) minimum micromolar concentration allowed
 %  param.externalNetFluxBounds:   ('original') =  
@@ -69,12 +67,15 @@
 
 %% processing for concentrations
 if ~isfield(param,'maxConc')
-    param.maxConc=1e4;
+    param.maxConc=inf;
 end
 if ~isfield(param,'minConc')
-    param.minConc=1e-4;
+    param.minConc=0;
 end
-
+% %assume units are in mMol
+% if ~isfield(param,'concUnit')
+%     param.concUnit = 10-3;
+% end
 
 if ~isfield(param,'externalNetFluxBounds')
     if isfield(model,'dcl') || isfield(model,'dcu')
@@ -87,8 +88,13 @@
 nMetabolitesPerRxn = sum(model.S~=0,1)';
 bool = nMetabolitesPerRxn>1 & ~model.SConsistentRxnBool;
 if any(bool)
-    warning('Exchange reactions involving more than one metabolite, check bounds on x - x0')
-    disp(model.rxns(bool))
+    warning([ int2str(nnz(bool)) ' stoichiometrically inconsistent reactions involving more than one metabolite, check bounds on x - x0'])
+    if nnz(bool)>10
+        ind=find(bool);
+        disp(model.rxns(ind(1:10)))
+    else
+        disp(model.rxns(bool))
+    end
 end
 
 if any(~model.SConsistentRxnBool)
@@ -253,10 +259,7 @@
     end
 end
 
-%assume concentrations are in uMol
-if ~isfield(model,'concUnit')
-    concUnit = 10-6;
-end
+
 
 % Define constants
 if isfield(model,'gasConstant') && isfield(model,'T')
@@ -303,12 +306,3 @@
         error('f must all be finite')
     end
 end
-
-% %lower and upper bounds on logarithmic concentration
-% pl = -log(param.maxConc*ones(m2,1));
-% if 1
-%     pu =  log(param.maxConc*ones(m2,1));
-% else
-%     %All potentials negative
-%     pu =  zeros(m2,1);
-% end
\ No newline at end of file
diff --git a/src/base/solvers/entropicFBA/processFluxConstraints.m b/src/base/solvers/entropicFBA/processFluxConstraints.m
index 65a61fb22d..3886c5d222 100644
--- a/src/base/solvers/entropicFBA/processFluxConstraints.m
+++ b/src/base/solvers/entropicFBA/processFluxConstraints.m
@@ -26,7 +26,7 @@
 %  param.printLevel:
 %  param.solver:    {'pdco',('mosek')}
 %  param.debug:     {(0),1} 1 = run in debug mode 
-%  param.method:    {('fluxes'),'fluxesConcentrations'} maximise entropy of fluxes (default) or also concentrations
+%  param.entropicFBAMethod:    {('fluxes'),'fluxesConcentrations'} maximise entropy of fluxes (default) or also concentrations
 %  param.maxUnidirectionalFlux: maximum unidirectional flux (1e5 by default)
 %  param.minUnidirectionalFlux: minimum unidirectional flux (zero by default)
 %  param.internalNetFluxBounds: ('original')   = use model.lb and model.ub to set the direction and magnitude of internal net flux bounds
@@ -56,7 +56,6 @@
 % Author(s): Ronan Fleming
 
 %% processing for fluxes
-
 if ~isfield(param,'maxUnidirectionalFlux')
     %try to set the maximum unidirectional flux based on the magnitude of the largest bound but dont have it greater than 1e5
     param.maxUnidirectionalFlux=min(1e5,max(abs(model.ub)));
@@ -74,8 +73,13 @@
     param.externalNetFluxBounds='original';
 end
 
-if ~isfield(param,'method')
-    param.method='fluxes';
+if ~isfield(param,'entropicFBAMethod')
+    if isfield(param,'method') && contains(param.method,'flux')
+        param.entropicFBAMethod=param.method;
+        param=rmfield(param,'method');
+    else
+        param.entropicFBAMethod='fluxes';
+    end
 end
 
 %find the maximal set of metabolites and reactions that are stoichiometrically consistent
@@ -98,8 +102,8 @@
     error('internalBounds replaced by other parameter options')
 end
 
-if isfield(param,'debug') && param.debug
-    solution_optimizeCbModel = optimizeCbModel(model);
+if isfield(param,'debug') && param.debug && 0
+    solution_optimizeCbModel = optimizeCbModel(model,param);
     switch solution_optimizeCbModel.stat
         case 0
             disp(solution_optimizeCbModel.origStat)
@@ -220,6 +224,7 @@
     vfl = model.vfl;
 else
     vfl = max(param.minUnidirectionalFlux,vl);
+    %vfl = ones(n,1)*param.minUnidirectionalFlux;
 end
 
 if any(vfl<0)
@@ -267,6 +272,7 @@
     vrl = model.vrl;
 else
     vrl = max(param.minUnidirectionalFlux,-vu);
+    %vrl = ones(n,1)*param.minUnidirectionalFlux;
 end
 if any(vrl<0)
     error('lower bound on reverse flux cannot be less than zero')
@@ -353,7 +359,7 @@
 end
 
 if ~isfield(model,'g') || isempty(model.g)
-    if isequal(param.method,'fluxes')
+    if isequal(param.entropicFBAMethod,'fluxes')
         model.g='one';
     else
         model.g='two';
diff --git a/src/base/solvers/entropicFBA/solveCobraEP.m b/src/base/solvers/entropicFBA/solveCobraEP.m
index 0f008254d1..9176859443 100644
--- a/src/base/solvers/entropicFBA/solveCobraEP.m
+++ b/src/base/solvers/entropicFBA/solveCobraEP.m
@@ -1,4 +1,4 @@
-function solution = solveCobraEP(EPproblem, varargin)
+function sol = solveCobraEP(EPproblem, varargin)
 % Solves the following optimisation problem:
 %
 % minimize   osense*(c.*d)'x + d.*x'(log(x) -1) + (1/2)*x'*Q*x
@@ -30,7 +30,7 @@
 %
 %
 % USAGE:
-%    solution = solveCobraEP(EPproblem, varargin)
+%    sol = solveCobraEP(EPproblem, varargin)
 %
 % INPUT:
 %    EPproblem:     Structure containing the following fields describing the EP problem to be solved
@@ -92,53 +92,57 @@
 %   optTol:         Optimality tolerance
 %
 % OUTPUT:
-%    solution:      Structure containing the following fields describing a LP solution:
-%                     * .obj:          Objective value
+%    sol:      Structure containing the following fields describing a LP sol:
+%                     *.obj:          Objective value
 %                     *.objLinear      osense*c'*x;
 %                     *.objEntropy     d.*x'*(log(x) -1);
 %                     *.objQuadratic   (1/2)*x'*Q*x;
-%                     * .full:         Primal solution vector
-%                     * .slack:        bl = A*x + s = bu
-%                     * .rcost:        Reduced costs, dual solution to :math:`lb <= x <= ub`
-%                     * .dual:         dual solution to constraints :math: `A*x ('E' | 'G' | 'L') b`
-%
-%                     * .solver:       Solver used to solve EP problem
+%                     *.v:  n+k ×1 double
+%                     *.vf: n × 1 double
+%                     *.vr: n × 1 double
+%                     *.vt: 1'*vt + 1'*vr
+%                     *.y_N: m x 1 double   dual sol to constraints :math: `A*x ('E' | 'G' | 'L') b`
+%                     *.z_dx: 0
+%                     *.z_vf: n × 1 double    dual sol to :math:`lb <= vr <= ub`
+%                     *.z_vr: n × 1 double    dual sol to :math:`lb <= vf <= ub`
+%                     *.z_vi: n × 1 double    dual sol to :math:`lb <= v <= ub`
+%                     *.z_v: n + k × 1 double dual sol to :math:`lb <= w <= ub`
 %                     * .stat:         Solver status in standardized form
 %                       * 0 - Infeasible problem
-%                       * 1 - Optimal solution
-%                       * 2 - Unbounded solution
-%                       * 3 - Almost optimal solution
+%                       * 1 - Optimal sol
+%                       * 2 - Unbounded sol
+%                       * 3 - Almost optimal sol
 %                       * -1 - Some other problem (timelimit, numerical problem etc)
 %                     * .origStat:         Original status returned by the specific solver
 %                     * .origStatText:     Original status text returned by the specific solver
 %                     * .time:         Solve time in seconds
-%
+%                     * .solver:       Solver used to solve EP problem
+%                     * .epmethod: solver method used e.g. 'CONIC'
+
 % OPTIONAL OUTPUT (from conic optimisation with mosek):
-%  solution.auxPrimal:  auxiliary primal variable
-%  solution.auxRcost:   dual to auxiliary primal variable
-%  solution.coneF:      affine constraint matrix
-%  solution.coneDual:   dual to affine constraints
-%  solution.dualNorm:   dual to the probability normalisation constraint
+%  sol.auxPrimal:  auxiliary primal variable
+%  sol.auxRcost:   dual to auxiliary primal variable
+%  sol.coneF:      affine constraint matrix
+%  sol.coneDual:   dual to affine constraints
+%  sol.dualNorm:   dual to the probability normalisation constraint
 %
 % OPTIONAL OUTPUT (from optimisation with pdco):
-%  solution.d1:  primal regularisation parameter, see pdco.m
-%  solution.d2:  dual regularisation parameter, see pdco.m
+%  sol.d1:  primal regularisation parameter, see pdco.m
+%  sol.d2:  dual regularisation parameter, see pdco.m
 %
 % EXAMPLE:
 %
-% NOTE: This code is a draft version released for the ELIXIR Fluxomic course and is not yet published and not to be redistributed without express permission of the author.
 %
-% Author(s): Ronan M.T. Fleming, 2021
+% Author(s): Ronan M.T. Fleming, 2024
 
-[problemTypeParams, solverOnlyParams] = parseSolverParameters('EP', varargin{:});
 
-if ~isfield(problemTypeParams,'debug')
-    problemTypeParams.debug = 1;
-end
+[problemTypeParams, solverParams] = parseSolverParameters('EP', varargin{:});
+param = mergeCobraParams(solverParams,problemTypeParams);
+clear problemTypeParams solverParams
 
-% Remove outer function specific parameters to avoid crashing solver interfaces
-% Default EP parameters are removed within solveCobraEP, so are not removed here
-solverOnlyParams = mosekParamStrip(solverOnlyParams);
+if ~isfield(param,'debug')
+    param.debug = 1;
+end
 
 if any(EPproblem.lb>EPproblem.ub)
     error('EPproblem.lb>EPproblem.ub');
@@ -153,7 +157,7 @@
     end
 end
 
-if isequal(problemTypeParams.solver,'mosek')
+if isequal(param.solver,'mosek')
     if ~(isfield(EPproblem,'blc') || isfield(EPproblem,'blc'))
         % blc <= A*x <= buc
         EPproblem.blc = EPproblem.b;
@@ -167,53 +171,38 @@
 end
 
 %% if in debug mode, test to see if the LP part of the problem is feasible
-if problemTypeParams.debug
-    switch problemTypeParams.solver
+if param.debug && ~isfield(param,'useTestVKSolution') && ~isfield(param,'VKproblem') %avoid LP if being called in driver_optimiseVKmodel
+    switch param.solver
         case 'pdco'
             solutionLP2 = solveCobraLP(EPproblem);
-            if problemTypeParams.printLevel>2
+            if param.printLevel>2
                 disp(solutionLP2)
             end
             
         case 'mosek'
-            %https://docs.mosek.com/8.1/toolbox/solving-linear.html
-            if ~isfield(problemTypeParams, 'MSK_DPAR_INTPNT_TOL_PFEAS')
-                solverOnlyParams.MSK_DPAR_INTPNT_TOL_PFEAS=problemTypeParams.feasTol;
-            end
-            if ~isfield(problemTypeParams, 'MSK_DPAR_INTPNT_TOL_DFEAS.')
-                solverOnlyParams.MSK_DPAR_INTPNT_TOL_DFEAS=problemTypeParams.feasTol;
-            end
+            [cmd,mosekParam] = setMosekParam(param);
 
-            %remove any fields with names that do not begin with 'MSK_'
-            solverOnlyParams = mosekParamStrip(solverOnlyParams);
-
-            [res] = msklpopt(EPproblem.c,EPproblem.A,EPproblem.blc,EPproblem.buc,EPproblem.lb,EPproblem.ub,solverOnlyParams,'minimize');
+            [res] = msklpopt(EPproblem.c,EPproblem.A,EPproblem.blc,EPproblem.buc,EPproblem.lb,EPproblem.ub,mosekParam,'minimize');
             
-            %If the feasibility tolerance is changed by the solverParams
-            %struct, this needs to be forwarded to the cobra Params for the
-            %final consistency test!
-            if isfield(problemTypeParams,'MSK_DPAR_INTPNT_TOL_PFEAS')
-                solverOnlyParams.feasTol = solverOnlyParams.MSK_DPAR_INTPNT_TOL_PFEAS;
-            end
-
             %parse mosek result structure
-            [solutionLP2.stat,solutionLP2.origStat,x,y,yl,yu,z,zl,zu,s,k,bas,pobjval,dobjval] = parseMskResult(res,EPproblem,solverOnlyParams,problemTypeParams.printLevel);
+           %[            stat,            origStat,x,y,yl,yu,z,zl,zu,k,basis,pobjval,dobjval] = parseMskResult(res)
+            [solutionLP2.stat,solutionLP2.origStat,x,y,yl,yu,z,zl,zu,k,  basis,pobjval,dobjval] = parseMskResult(res);
         
             switch solutionLP2.stat
                 case 0
-                    solution = solutionLP2;
-                    message = ['solveCobraEP: LP part of EPproblem is infeasible according to solveCobraLP with ' problemTypeParams.solver '.'];
+                    sol = solutionLP2;
+                    message = ['solveCobraEP: LP part of EPproblem is infeasible according to solveCobraLP with ' param.solver '.'];
                     warning(message)
                     
                     return
                 case 2
-                    solution = solutionLP2;
-                    message = ['solveCobraEP: LP part of EPproblem is unbounded according to solveCobraLP with ' problemTypeParams.solver '.'];
+                    sol = solutionLP2;
+                    message = ['solveCobraEP: LP part of EPproblem is unbounded according to solveCobraLP with ' param.solver '.'];
                     warning(message)
                     
                     return
                 case 1
-                    message =['solveCobraEP: LP part of EPproblem is feasible according to solveCobraLP with ' problemTypeParams.solver '.'];
+                    message =['solveCobraEP: LP part of EPproblem is feasible according to solveCobraLP with ' param.solver '.'];
                     fprintf('%s\n',message)
                 otherwise
                     error('inconclusive solveCobraLP')
@@ -277,7 +266,7 @@
 
 [mlt,nlt]=size(A);
 
-switch problemTypeParams.solver
+switch param.solver
     case 'pdco'
         % solves optimization problems of the form
         %
@@ -334,38 +323,38 @@
             end
         end
                 
-        if isfield(solverOnlyParams,'d1')
-            d1 = solverOnlyParams.d1;
+        if isfield(param,'d1')
+            d1 = param.d1;
         else
             d1 = 1e-4;
         end
-        if isfield(solverOnlyParams,'d2')
-            d2 = solverOnlyParams.d2;
+        if isfield(param,'d2')
+            d2 = param.d2;
         else
             d2 = 1e-4;
         end
-        if isfield(solverOnlyParams,'x0')
-            x0 = solverOnlyParams.x0;
+        if isfield(param,'x0')
+            x0 = param.x0;
         else
             x0 = ones(size(Aeq,2),1);
         end
-        if isfield(solverOnlyParams,'y0')
-            y0 = solverOnlyParams.y0;
+        if isfield(param,'y0')
+            y0 = param.y0;
         else
             y0 = ones(size(Aeq,1),1);
         end
-        if isfield(solverOnlyParams,'z0')
-            z0 = solverOnlyParams.z0;
+        if isfield(param,'z0')
+            z0 = param.z0;
         else
             z0 = ones(size(Aeq,2),1);
         end
-        if isfield(solverOnlyParams,'xsize')
-            xsize = solverOnlyParams.xsize;
+        if isfield(param,'xsize')
+            xsize = param.xsize;
         else
             xsize = 1;
         end
-        if isfield(solverOnlyParams,'zsize')
-            zsize = solverOnlyParams.zsize;
+        if isfield(param,'zsize')
+            zsize = param.zsize;
         else
             zsize = 1;
         end
@@ -374,8 +363,8 @@
         options = pdcoSet;
         %options.mu0       = 1; %very small only for entropy function
         options.mu0       = 0; %pdco chooses its own
-        options.FeaTol    = problemTypeParams.feasTol;
-        options.OptTol    = problemTypeParams.optTol;
+        options.FeaTol    = param.feasTol;
+        options.OptTol    = param.optTol;
         %   If getting linesearch failures, slacken tolerances
         %   i.e. Linesearch failed (nf too big)
         %options.FeaTol    = 1e-6; %%Ecoli core working at 1e-7
@@ -443,7 +432,7 @@
         else
             objHandle = @(x) entropyObj(x,ceq,deq);
         end
-        options.Print = problemTypeParams.printLevel-1;
+        options.Print = param.printLevel-1;
         
         saveAndDebug=0;
         if saveAndDebug
@@ -457,7 +446,7 @@
         logx(deq~=0) = reallog(x(deq~=0));     % error if negative
         grad = ceq + deq.*logx;
         
-        if problemTypeParams.printLevel > 2 || problemTypeParams.debug
+        if param.printLevel > 2 || param.debug
             % determine the residuals
             fprintf('\n%s\n','KKT with pdco signs:')
             fprintf('%8.2g %s\n',norm(Aeq*x - beq,inf), '|| Aeq*x - beq ||_inf');
@@ -473,15 +462,16 @@
             return
         end
         
-        solution.time = toc;
-        
-        % inform = 0 if a solution is found;
+        sol.time = toc;
+        sol.epmethod = options.Method;
+
+        % inform = 0 if a sol is found;
         %        = 1 if too many iterations were required;
         %        = 2 if the linesearch failed too often;
         %        = 3 if the step lengths became too small;
         %        = 4 if Cholesky said ADDA was not positive definite.
         if (inform == 0)
-            solution.stat = 1;
+            sol.stat = 1;
             
 
             if ~any(csense == 'L' | csense == 'G')
@@ -495,71 +485,73 @@
                 y(csense == 'G') = -y(csense == 'G');
             end
             
-            solution.slack = slack;
-            solution.full = x(1:nlt,1);
-            solution.dual = -y;
-            solution.rcost = -z(1:nlt,1);
-            solution.origStat = inform;
+            sol.slack = slack;
+            sol.full = x(1:nlt,1);
+            sol.dual = -y;
+            sol.rcost = -z(1:nlt,1);
+            sol.origStat = inform;
 
             %objective
             logx = zeros(size(A,2),1);
-            logx(d~=0) = reallog(solution.full(d~=0));     % error if negative
+            logx(d~=0) = reallog(sol.full(d~=0));     % error if negative
             if isfield(EPproblem,'Q')
-                solution.obj  = c'*solution.full + (d.*solution.full)'*logx + (1/2)*solution.full'*EPproblem.Q*solution.full;
-                grad = c + d.*logx + EPproblem.Q*solution.full;
+                sol.obj  = c'*sol.full + (d.*sol.full)'*logx + (1/2)*sol.full'*EPproblem.Q*sol.full;
+                grad = c + d.*logx + EPproblem.Q*sol.full;
             else
-                solution.obj  = c'*solution.full + (d.*solution.full)'*logx;
+                sol.obj  = c'*sol.full + (d.*sol.full)'*logx;
                 grad = c + d.*logx;
             end
             Aty = -A'*y;
             
-            if problemTypeParams.printLevel > 2 || problemTypeParams.debug
+            if param.printLevel > 2 || param.debug
                 fprintf('\n%s\n','KKT with Rockafellar signs:')
-                fprintf('%8.2g %s\n',norm(A*solution.full + solution.slack - b,inf), '|| A*x + s - b ||_inf');
-                fprintf('%8.2g %s\n',norm(A*solution.full + solution.slack - b - (d2^2)*solution.dual,inf), '|| A*x + s - b - (d2^2)*y ||_inf');
-                res2 = grad  + A'*solution.dual + solution.rcost;
+                fprintf('%8.2g %s\n',norm(A*sol.full + sol.slack - b,inf), '|| A*x + s - b ||_inf');
+                fprintf('%8.2g %s\n',norm(A*sol.full + sol.slack - b - (d2^2)*sol.dual,inf), '|| A*x + s - b - (d2^2)*y ||_inf');
+                res2 = grad  + A'*sol.dual + sol.rcost;
                 fprintf('%8.2g %s\n',norm(res2,inf), '|| grad + A''*y + z ||_inf');
                 
-                fprintf('%8.2g %s\n',norm(grad  + A'*solution.dual + solution.rcost - (d1^2)*solution.full,inf), '|| grad + A''*y + z - (d1^2)*x ||_inf');
-                if problemTypeParams.debug
-                    %res2 = grad  + Aty + solution.rcost;
-                    res2 = grad  + A'*solution.dual + solution.rcost;
-                    solution.T = table(res2,c,d.*logx, -A'*y, solution.rcost,...
+                fprintf('%8.2g %s\n',norm(grad  + A'*sol.dual + sol.rcost - (d1^2)*sol.full,inf), '|| grad + A''*y + z - (d1^2)*x ||_inf');
+                if param.debug
+                    %res2 = grad  + Aty + sol.rcost;
+                    res2 = grad  + A'*sol.dual + sol.rcost;
+                    sol.T = table(res2,c,d.*logx, -A'*y, sol.rcost,...
                         'VariableNames',{'total','c','dlogx','Aty','z'});
                 end
                 if any(~isfinite(res2))
                     warning('Infinite variables in dual optimality condition')
-                    solution.Tinf=solution.T(~isfinite(solution.T.total),:);
-                    ind = find(~isfinite(solution.T.total));
-                    solution.Tinf.ind = ind;
+                    sol.Tinf=sol.T(~isfinite(sol.T.total),:);
+                    ind = find(~isfinite(sol.T.total));
+                    sol.Tinf.ind = ind;
                 end
             end
             
-            solution.d1=d1;
-            solution.d2=d2;
+            sol.d1=d1;
+            sol.d2=d2;
             
         elseif (inform == 1 || inform == 2 || inform == 3)
-            solution.stat = 0;
-            solution.obj = NaN;
+            sol.stat = 0;
+            sol.obj = NaN;
         else
-            solution.stat = -1;
-            solution.obj = NaN;
+            sol.stat = -1;
+            sol.obj = NaN;
         end
-        solution.origStat = inform;
+        sol.origStat = inform;
         
         %update parameters for testing optimality criterion
-        problemTypeParams.feasTol = options.FeaTol;
-        problemTypeParams.optTol = options.OptTol;
-        %                     * .full:         Full LP solution vector
+        param.feasTol = options.FeaTol;
+        param.optTol = options.OptTol;
+        %                     * .full:         Full LP sol vector
         %                     * .obj:          Objective value
-        %                     * .rcost:        Reduced costs, dual solution to :math:`lb <= v <= ub`
-        %                     * .dual:         dual solution to `A*v ('E' | 'G' | 'L') b`
+        %                     * .rcost:        Reduced costs, dual sol to :math:`lb <= v <= ub`
+        %                     * .dual:         dual sol to `A*v ('E' | 'G' | 'L') b`
         %                     * .solver:       Solver used to solve LP problem
         %                     * .algorithm:    Algorithm used by solver to solve LP problem
         %                     * .stat:         Solver status in standardized form
                
     case 'mosek'
         %%
+        [cmd, mosekParam] = setMosekParam(param);
+
         %         https://docs.mosek.com/modeling-cookbook/expo.html
         %         https://docs.mosek.com/modeling-cookbook/qcqo.html#conic-reformulation
         %         min  (d.*x)'*(log(x./y) + c)  + (1/2)*x'*Q*x
@@ -669,8 +661,8 @@
 
         % Specify conic part of the problem
         % https://docs.mosek.com/9.2/toolbox/data-types.html#cones
-        if problemTypeParams.printLevel>1 || problemTypeParams.debug
-            [~, res] = mosekopt('symbcon');
+        if param.printLevel>1 || param.debug
+            [~, res] = mosekopt('symbcon echo(0)');
         else
             [~, res] = mosekopt('symbcon echo(0)');
         end
@@ -809,75 +801,36 @@
             mosekopt('write(problem.opf)',prob)
         end
         
-        %set default mosek parameters for this type of problem
-        paramMosek=mosekParamSetEFBA;
-        
-        if ~isfield(solverOnlyParams,'MSK_DPAR_INTPNT_CO_TOL_PFEAS')
-            if isfield(solverOnlyParams,'MSK_DPAR_INTPNT_CO_TOL_PFEAS')
-                paramMosek.MSK_DPAR_INTPNT_CO_TOL_PFEAS = solverOnlyParams.feasTol;
-            else
-                paramMosek.MSK_DPAR_INTPNT_CO_TOL_PFEAS = problemTypeParams.feasTol;
-            end
-        end
-        if ~isfield(solverOnlyParams,'MSK_DPAR_INTPNT_CO_TOL_DFEAS')
-            if isfield(solverOnlyParams,'MSK_DPAR_INTPNT_CO_TOL_DFEAS')
-                paramMosek.MSK_DPAR_INTPNT_CO_TOL_DFEAS = solverOnlyParams.optTol;
-            else
-                paramMosek.MSK_DPAR_INTPNT_CO_TOL_DFEAS = problemTypeParams.optTol;
-            end
-        end
-        
-        % only set the print level if not already set via solverParams structure
-        if ~isfield(solverOnlyParams, 'MSK_IPAR_LOG')
-            switch problemTypeParams.printLevel
-                case 0
-                    echolev = 0;
-                case 1
-                    echolev = 3;
-                case 2
-                    paramMosek.MSK_IPAR_LOG_INTPNT = 1;
-                    paramMosek.MSK_IPAR_LOG_SIM = 1;
-                    echolev = 3;
-                otherwise
-                    echolev = 0;
-            end
-            if echolev == 0
-                paramMosek.MSK_IPAR_LOG = 0;
-                cmd = ['minimize echo(' int2str(echolev) ')'];
-            else
-                cmd = 'minimize';
-            end
-            
-        end
-        %overide if in debug mode
-        if problemTypeParams.debug
-            cmd = 'minimize';
-        end
-            
-        if problemTypeParams.debug && 0
-            probBeforeMosekopt = prob;
-            save('probBeforeMosekopt','probBeforeMosekopt');
+
+        if isfield(param,'saveProb') && param.saveProb
+            formattedTime = datestr(now, 'yyyymmddHHMMSS');
+            EP.cmd=cmd;
+            EP.prob=prob;
+            EP.param=mosekParam;
+            save([formattedTime '_EP_probBeforeMosekopt'],"EP");
         end
 
-        %param = updateStructData(param,solverParams);
-        
         %call mosek exponential cone solver
         tic;
         if 0
             %default
             [~,res]=mosekopt('minimize',prob);
         else
-            [~,res]=mosekopt(cmd,prob,paramMosek);
+            [~,res]=mosekopt(cmd,prob,mosekParam);
         end
-        solution.time = toc;
-        
+        sol.time = toc;
+        if isfield(param,'epmethod')
+            sol.epmethod = param.epmethod;
+        else
+            sol.epmethod = 'FREE';
+        end
+
         %parse mosek result structure      
-       %[stat,origStat,x,y,yl,yu,z,zl,zu,k,basis,pobjval,dobjval] = parseMskResult(res,solverParams,printLevel)
-        [stat,origStat,x,y,yl,yu,z,zl,zu,k,bas,pobjval,dobjval] = parseMskResult(res,solverOnlyParams,problemTypeParams.printLevel);
-       %[stat,origStat,x,y,yl,yu,z,zl,zu,k,basis,pobjval,dobjval] = parseMskResult(res,solverOnlyParams,printLevel)
+        [stat,origStat,x,y,yl,yu,z,zl,zu,k,basis,pobjval,dobjval] = parseMskResult(res);
 
-        solution.stat = stat;
-        solution.origStat = origStat;
+        sol.stat = stat;
+        sol.origStat = origStat;
+        sol.obj = pobjval;
         switch stat
             case 1
                 %check for zeros in variables within entropy functions
@@ -895,20 +848,20 @@
                 sbl = prob.a*x - prob.blc;
                 sbu = prob.buc - prob.a*x;
                 s = sbu - sbl; %TODO -double check this
-                if problemTypeParams.printLevel>1
+                if param.printLevel>1
                     fprintf('%8.2g %s\n',min(sbl), ' min(sbl) = min(A*x - bl), (should be positive)');
                     fprintf('%8.2g %s\n',min(sbu), ' min(sbu) = min(bu - A*x), (should be positive)');
                 end
                 
-                if problemTypeParams.printLevel > 1
+                if param.printLevel > 1
                     % Problem definition here: https://docs.mosek.com/9.2/toolbox/prob-def-affine-conic.html
                     fprintf('%s\n','Optimality conditions (numerical)')
-                    % Guide to interpreting the solution summary: https://docs.mosek.com/9.2/toolbox/debugging-log.html#continuous-problem
+                    % Guide to interpreting the sol summary: https://docs.mosek.com/9.2/toolbox/debugging-log.html#continuous-problem
                     fprintf('%8.2g %s\n',norm(prob.a(prob.blc==prob.buc,:)*x - prob.blc(prob.blc==prob.buc),inf), '|| A*x - b ||_inf');
                     val = norm(prob.c - prob.a'*y - z - prob.f'*k,inf);
                     fprintf('%8.2g %s\n',val, '|| c - A''*y - z - F''*k ||_inf');
-                    if val>1e-6 || problemTypeParams.debug
-                        solution.T0 = table(prob.c - prob.a'*y - z - prob.f'*k,prob.c, prob.a'*y, z,prob.f'*k,'VariableNames',{'tot','c','Aty','z','Ftdoty'});
+                    if val>1e-6 || param.debug
+                        sol.T0 = table(prob.c - prob.a'*y - z - prob.f'*k,prob.c, prob.a'*y, z,prob.f'*k,'VariableNames',{'tot','c','Aty','z','Ftdoty'});
                     end
                     %fprintf('%8.2g %s\n',norm(prob.c - prob.f'*s,inf), '|| c - F''k ||_inf');
                     fprintf('%8.2g %s\n',norm(-y + res.sol.itr.slc - res.sol.itr.suc,inf), '|| -y + res.sol.itr.slc - res.sol.itr.suc ||_inf');
@@ -926,14 +879,14 @@
                 
                 %check with the original order of the affine cone constraints
                 val = norm(prob.c - prob.a'*y - z - F'*y_K,inf);
-                if problemTypeParams.printLevel > 1
+                if param.printLevel > 1
                     fprintf('%8.2g %s\n',val, '|| c - A''*y - z - F''*y_K ||_inf');
                 end
-                if val>1e-6 || problemTypeParams.debug
-                    solution.T = table(prob.c - prob.a'*y - z - F'*y_K,prob.c, prob.a'*y, z,prob.f'*k,F'*y_K,'VariableNames',{'tot','c','Aty','z','Ftdoty','Fty_K'});
+                if val>1e-6 || param.debug
+                    sol.T = table(prob.c - prob.a'*y - z - F'*y_K,prob.c, prob.a'*y, z,prob.f'*k,F'*y_K,'VariableNames',{'tot','c','Aty','z','Ftdoty','Fty_K'});
                 end
                 
-                if problemTypeParams.printLevel > 1
+                if param.printLevel > 1
                     x1 = F(1:nCone,:)*x;
                     x2 = F(nCone+1:2*nCone,:)*x;
                     x3 = F(2*nCone+1:3*nCone,:)*x;
@@ -973,40 +926,40 @@
                     fprintf('%7.2g\t%s\n',min(y1_K(1:nExpCone) + y3_K(1:nExpCone).*exp(y2_K(1:nExpCone)./y3_K(1:nExpCone))/exp(1)), 'min(y1_k + y3_k.*exp(y2_K./y3_K)/exp(1))  >= 0');
                 end
                 
-                solution.full = x(1:size(A,2));
+                sol.full = x(1:size(A,2));
                 %switch to Rockafellar signs
-                solution.dual = -y(1:size(A,1));
-                solution.dualNorm = -y(size(A,1)+1:size(A,1)+p);
-                solution.rcost = -z(1:size(A,2)+p);
-                solution.slack = s;
+                sol.dual = -y(1:size(A,1));
+                sol.dualNorm = -y(size(A,1)+1:size(A,1)+p);
+                sol.rcost = -z(1:size(A,2)+p);
+                sol.slack = s;
                 
                 %need to zero out the NaN due to log(0) for some variables
-                logSolutionFull = real(log(solution.full));
+                logSolutionFull = real(log(sol.full));
                 logSolutionFull(~isfinite(logSolutionFull))=0;
                 if isfield(EPproblem,'Q')
-                    solution.obj = EPproblem.c'*solution.full + (EPproblem.d.*solution.full)'*(logSolutionFull -1) + (1/2)*solution.full'*EPproblem.Q*solution.full;
-                    solution.objLinear = EPproblem.c'*solution.full;
-                    solution.objEntropy = -(EPproblem.d.*solution.full)'*(logSolutionFull -1);
-                    solution.objQuadratic = (1/2)*solution.full'*EPproblem.Q*solution.full;
+                    sol.obj = EPproblem.c'*sol.full + (EPproblem.d.*sol.full)'*(logSolutionFull -1) + (1/2)*sol.full'*EPproblem.Q*sol.full;
+                    sol.objLinear = EPproblem.c'*sol.full;
+                    sol.objEntropy = -(EPproblem.d.*sol.full)'*(logSolutionFull -1);
+                    sol.objQuadratic = (1/2)*sol.full'*EPproblem.Q*sol.full;
                 else
-                    solution.obj = EPproblem.c'*solution.full + (EPproblem.d.*solution.full)'*(logSolutionFull -1);
-                    solution.objLinear = EPproblem.c'*solution.full;
-                    solution.objEntropy = -(EPproblem.d.*solution.full)'*(logSolutionFull -1);
-                    solution.objQuadratic = 0;
+                    sol.obj = EPproblem.c'*sol.full + (EPproblem.d.*sol.full)'*(logSolutionFull -1);
+                    sol.objLinear = EPproblem.c'*sol.full;
+                    sol.objEntropy = -(EPproblem.d.*sol.full)'*(logSolutionFull -1);
+                    sol.objQuadratic = 0;
 
                 end
                 
-                posRcost = solution.rcost>0;
-                negRcost = solution.rcost<0;
+                posRcost = sol.rcost>0;
+                negRcost = sol.rcost<0;
                 blx = prob.blx(1:size(A,2));
                 bux = prob.bux(1:size(A,2));
-                solution.lagRcost = sum(solution.rcost(negRcost)'*blx(negRcost) + solution.rcost(posRcost)'*bux(posRcost));
+                sol.lagRcost = sum(sol.rcost(negRcost)'*blx(negRcost) + sol.rcost(posRcost)'*bux(posRcost));
                 
                 %pass back the F matrix to check biochemical optimality criteria
-                solution.coneF = F;
-                solution.auxPrimal = x(size(A,2)+p+1:end);
-                solution.auxRcost = -z(size(A,2)+p+1:end);
-                solution.coneDual = -y_K;
+                sol.coneF = F;
+                sol.auxPrimal = x(size(A,2)+p+1:end);
+                sol.auxRcost = -z(size(A,2)+p+1:end);
+                sol.coneDual = -y_K;
                 
                 % variable to determine the residual 1
                 b = prob.blc;
@@ -1021,146 +974,111 @@
                 k = NaN*ones(size(prob.f,1),1);
         end              
     otherwise
-        error([problemTypeParams.solver ' is an unrecognised solver'])
+        error([param.solver ' is an unrecognised solver'])
 end
 
-switch solution.stat
+sol.solver=param.solver;
+
+switch sol.stat
     case 0
-        switch problemTypeParams.solver
+        switch param.solver
             case 'pdco'
                 %infeasible, debug the situtation
-                disp(solution.origStat)
-                %solution.origStat: 'PRIMAL_INFEASIBLE_CER'
+                disp(sol.origStat)
+                %sol.origStat: 'PRIMAL_INFEASIBLE_CER'
                 solutionLP = solveCobraLP(EPproblem);
                 statLP=solutionLP.stat;
                 
             case 'mosek'
-                %https://docs.mosek.com/8.1/toolbox/solving-linear.html
-                if ~isfield(problemTypeParams, 'MSK_DPAR_INTPNT_TOL_PFEAS')
-                    solverOnlyParams.MSK_DPAR_INTPNT_TOL_PFEAS=problemTypeParams.feasTol;
-                end
-                if ~isfield(problemTypeParams, 'MSK_DPAR_INTPNT_TOL_DFEAS.')
-                    solverOnlyParams.MSK_DPAR_INTPNT_TOL_DFEAS=problemTypeParams.feasTol;
-                end
-                %If the feasibility tolerance is changed by the solverParams
-                %struct, this needs to be forwarded to the cobra Params for the
-                %final consistency test!
-                if isfield(problemTypeParams,'MSK_DPAR_INTPNT_TOL_PFEAS')
-                    solverOnlyParams.feasTol = solverOnlyParams.MSK_DPAR_INTPNT_TOL_PFEAS;
-                end
+                [cmd,mosekParam] = setMosekParam(param);
 
-                
-                % only set the print level if not already set via solverParams structure
-                if ~isfield(solverOnlyParams, 'MSK_IPAR_LOG')
-                    switch problemTypeParams.printLevel
-                        case 0
-                            echolev = 0;
-                        case 1
-                            echolev = 3;
-                        case 2
-                            solverOnlyParams.MSK_IPAR_LOG_INTPNT = 1;
-                            solverOnlyParams.MSK_IPAR_LOG_SIM = 1;
-                            echolev = 3;
-                        otherwise
-                            echolev = 0;
-                    end
-                end
-                if echolev == 0
-                    solverOnlyParams.MSK_IPAR_LOG = 0;
-                    cmd = ['minimize echo(' int2str(echolev) ')'];
-                else
-                    cmd = 'minimize';
-                end
-                [res] = msklpopt(EPproblem.c,EPproblem.A,EPproblem.blc,EPproblem.buc,EPproblem.lb,EPproblem.ub,solverOnlyParams,cmd);
-                
-                %[stat,origStat,x,y,yl,yu,z,zl,zu,s,k,bas,pobjval,dobjval] = parseMskResult(res,prob,problemTypeParams.printLevel);
-                [statLP,origStat,x,y,yl,yu,z,zl,zu,k,bas,pobjval,dobjval] = parseMskResult(res,solverOnlyParams,problemTypeParams.printLevel);
-                
-   
+                [res] = msklpopt(EPproblem.c,EPproblem.A,EPproblem.blc,EPproblem.buc,EPproblem.lb,EPproblem.ub,mosekParam,cmd);
+
+                [statLP,origStat,x,y,yl,yu,z,zl,zu,k,basis,pobjval,dobjval] = parseMskResult(res);
         end
         
         switch statLP
             case 1
-                message =['solveCobraEP: EPproblem with ' problemTypeParams.solver ' is infeasible, but corresponding LPproblem is feasible according to solveCobraLP with ' problemTypeParams.solver];
+                message =['solveCobraEP: EPproblem with ' param.solver ' is infeasible, but corresponding LPproblem is feasible according to solveCobraLP with ' param.solver];
                 warning(message)
             otherwise
-                message = ['solveCobraEP: EPproblem with ' problemTypeParams.solver ' is infeasible, because corresponding LPproblem is infeasible according to solveCobraLP with ' problemTypeParams.solver];
+                message = ['solveCobraEP: EPproblem with ' param.solver ' is infeasible, because corresponding LPproblem is infeasible according to solveCobraLP with ' param.solver];
                 warning(message)
         end
         if exist('messages','var')
-            if isfield(solution,'messages')
-                solution.messages = [messages;solution.messages;message];
+            if isfield(sol,'messages')
+                sol.messages = [messages;sol.messages;message];
             else
-                solution.messages = [messages;message];
+                sol.messages = [messages;message];
             end
         else
-            solution.messages = cellstr(message);
+            sol.messages = cellstr(message);
         end
       
     case 1
         % check the optimality conditions for various solvers
-        if ~isempty(solution.slack) && ~isempty(solution.full)
+        if ~isempty(sol.slack) && ~isempty(sol.full)
             % determine the residual 1
-            switch problemTypeParams.solver
+            switch param.solver
                 case 'pdco'
                     feasTol = 1e-3;
-                    res1 = A*solution.full + solution.slack - b;
+                    res1 = A*sol.full + sol.slack - b;
                     res1(~isfinite(res1))=0;
                 case 'mosek'
-                    feasTol = problemTypeParams.feasTol * 1e2;
-                    res1 = A(blc==buc,:)*solution.full - blc(blc==buc);
+                    feasTol = param.feasTol * 1e2;
+                    res1 = A(blc==buc,:)*sol.full - blc(blc==buc);
             end
             tmp1 = norm(res1, inf);
             
             % evaluate the optimality condition 1
             if tmp1 > feasTol
-                if strcmp(problemTypeParams.solver,'pdco')
-                    res1b = norm(A*solution.full + solution.slack - b + (d2^2)*y,inf);
+                if strcmp(param.solver,'pdco')
+                    res1b = norm(A*sol.full + sol.slack - b + (d2^2)*y,inf);
                     tmp1b = norm(res1b, inf);
                     if tmp1b > feasTol
                         displayError = 1;
                     else
                         displayError = 0;
-                        warning(['[' problemTypeParams.solver '] Primal optimality condition in solveCobraEP only approximately satisfied, residual = ' num2str(tmp1) ', regularised residual = ' num2str(tmp1b) ', while problem feasTol = ' num2str(feasTol) '.  origStat = ' solution.origStat])
+                        warning(['[' param.solver '] Primal optimality condition in solveCobraEP only approximately satisfied, residual = ' num2str(tmp1) ', regularised residual = ' num2str(tmp1b) ', while problem feasTol = ' num2str(feasTol) '.  origStat = ' sol.origStat])
                     end
                 else
                     %TODO - debug why solver reporting optimal but unscaled seems less so.
                     displayError = 0;
                 end
                 if displayError
-                    %disp(solution.origStat)
-                    fprintf('%s\n',['[' problemTypeParams.solver '] Primal optimality condition in solveCobraEP not satisfied, residual = ' num2str(tmp1) ', while problem feasTol = ' num2str(feasTol) '.  origStat = ' solution.origStat])
+                    %disp(sol.origStat)
+                    fprintf('%s\n',['[' param.solver '] Primal optimality condition in solveCobraEP not satisfied, residual = ' num2str(tmp1) ', while problem feasTol = ' num2str(feasTol) '.  origStat = ' sol.origStat])
                 end
             else
-                if problemTypeParams.printLevel > 0
-                    fprintf(['\n > [' problemTypeParams.solver '] Primal optimality condition in solveCobraEP satisfied.']);
+                if param.printLevel > 0
+                    fprintf(['\n > [' param.solver '] Primal optimality condition in solveCobraEP satisfied.']);
                 end
             end
         end
         
         %gradient may differ depending on the solver
-        res2 = grad  + Aty + solution.rcost;
+        res2 = grad  + Aty + sol.rcost;
         tmp2 = norm(res2, inf);
         
         if 0
-            optTol = problemTypeParams.optTol * 1e2;
+            optTol = param.optTol * 1e2;
         else
             optTol = 5e-5;
         end
         % evaluate the optimality condition 2
         if tmp2 > optTol
-            disp(solution.origStat)
-            if ~(length(A)==1 && strcmp(problemTypeParams.solver,'pdco')) %todo, why does pdco choke on small A?
-                warning(['[' problemTypeParams.solver '] Dual   optimality condition in solveCobraEP not satisfied, residual = ' num2str(tmp2) ', while problem optTol = ' num2str(optTol)])
+            disp(sol.origStat)
+            if ~(length(A)==1 && strcmp(param.solver,'pdco')) %todo, why does pdco choke on small A?
+                warning(['[' param.solver '] Dual   optimality condition in solveCobraEP not satisfied, residual = ' num2str(tmp2) ', while problem optTol = ' num2str(optTol)])
             end
         else
-            if problemTypeParams.printLevel > 0
-                fprintf(['\n > [' problemTypeParams.solver '] Dual   optimality condition in solveCobraEP satisfied.\n']);
+            if param.printLevel > 0
+                fprintf(['\n > [' param.solver '] Dual   optimality condition in solveCobraEP satisfied.\n']);
             end
         end
 end
 
-solution.solver=problemTypeParams.solver;
+
 
 end
 
diff --git a/src/base/solvers/getSetSolver/changeCobraSolver.m b/src/base/solvers/getSetSolver/changeCobraSolver.m
index ccc5c7839b..24e4863f4c 100644
--- a/src/base/solvers/getSetSolver/changeCobraSolver.m
+++ b/src/base/solvers/getSetSolver/changeCobraSolver.m
@@ -1,4 +1,4 @@
-function [solverOK, solverInstalled] = changeCobraSolver(solverName, solverType, printLevel, validationLevel)
+function [solverOK, solverInstalled] = changeCobraSolver(solverName, problemType, printLevel, validationLevel)
 % Changes the Cobra Toolbox optimization solver(s)
 %
 % USAGE:
@@ -7,10 +7,10 @@
 %
 % INPUTS:
 %    solverName:           Solver name
-%    solverType:           Solver type ('LP' by default)
+%    problemType:          Problem type ('LP' by default)
 %                          (a) One of the following: `LP` `MILP`, `QP`, `MIQP` 'EP', 'CLP'
 %                          (b) 'all' attempts to change all applicable solvers to solverName.  This is purely a shorthand convenience.
-%                          (c) Cell array of solverTypes, e.g. {'LP','QP'}   
+%                          (c) Cell array of problemTypes, e.g. {'LP','QP'}   
 %    printLevel:           verbose level
 %
 %                           *   if `0`, warnings and errors are silenced
@@ -21,7 +21,7 @@
 %
 %                           *   `-1`: assign only the global variable. Do not assign any path.
 %                           *   `0`: adjust solver paths but don't validate the solver
-%                           *   `1`: validate but remove outputs  (default)
+%                           *   `1`: validate but remove outputs, silent  (default)
 %                           *   `2`: validate and keep any outputs
 %
 % OUTPUT:
@@ -207,7 +207,7 @@
 solverInstalled = true;
 
 if validationLevel == -1
-    switch solverType
+    switch problemType
         case 'LP'
             CBT_LP_SOLVER = solverName;
         case 'QP'
@@ -299,9 +299,9 @@
 end
 
 if nargin < 2
-    solverType = 'LP';
+    problemType = 'LP';
 else
-    solverType = upper(solverType);
+    problemType = upper(problemType);
 end
 
 % print an error message if the solver is not supported
@@ -333,12 +333,12 @@
 end
 
 % Attempt to set the user provided solver for all optimization problem types
-if strcmpi(solverType, 'all')
+if strcmpi(problemType, 'all')
     solvedProblems = SOLVERS.(solverName).type;
     for i = 1:length(solvedProblems)
         %this looks to be self referential (calling changeCobraSolver in
         %changeCobraSolver)
-        [solverOK,solverInstalled] = changeCobraSolver(solverName, solvedProblems{i}, printLevel);
+        [solverOK,solverInstalled] = changeCobraSolver(solverName, solvedProblems{i}, printLevel,validationLevel);
         if printLevel > 0
             fprintf([' > changeCobraSolver: Solver for ', solvedProblems{i}, ' problems has been set to ', solverName, '.\n']);
         end
@@ -360,24 +360,24 @@
 
 % check if the given solver is able to solve the given problem type.
 solverOK = false;
-if ~contains(solverType, OPT_PROB_TYPES)
+if ~contains(problemType, OPT_PROB_TYPES)
     %This is not done during init, so at this point, the solver is already
     %checked for installation
     solverInstalled = SOLVERS.(solverName).installed;
     if printLevel > 0
-        error('changeCobraSolver: %s  problems cannot be solved in The COBRA Toolbox', solverType);
+        error('changeCobraSolver: %s  problems cannot be solved in The COBRA Toolbox', problemType);
     else
         return
     end
 end
 
 % check if the given solver is able to solve the given problem type.
-if ~contains(solverType, SOLVERS.(solverName).type)
+if ~contains(problemType, SOLVERS.(solverName).type)
     %This is not done during init, so at this point, the solver is already
     %checked for installation
     solverInstalled = SOLVERS.(solverName).installed;
     if printLevel > 0
-        error('Solver %s cannot solve %s problems', solverName, solverType);
+        error('Solver %s cannot solve %s problems', solverName, problemType);
     else
         return
     end
@@ -528,13 +528,13 @@
     if validationLevel > 0
         cwarn = warning;
         warning('off');
-        eval(['oldval = CBT_', solverType, '_SOLVER;']);
-        eval(['CBT_', solverType, '_SOLVER = solverName;']);
+        eval(['oldval = CBT_', problemType, '_SOLVER;']);
+        eval(['CBT_', problemType, '_SOLVER = solverName;']);
         % validate with a simple problem.
-        if strcmp(solverName,'mosek') && strcmp(solverType,'CLP') || strcmp(solverType,'all')
-            problem = struct('A',[0 1],'b',0,'c',[1;1],'osense',-1,'lb',[0;0],'ub',[0;0],'csense','E','vartype',['C';'I'],'x0',[0;0]);  
+        if strcmp(solverName,'mosek') && strcmp(problemType,'CLP') || strcmp(problemType,'all')
+            problem = struct('A',[0 1],'b',0,'c',[1;1],'osense',-1,'lb',[0;0],'ub',[0;0],'csense','E','vartype',['C';'I'],'x0',[0;0],'names',[]);  
         else
-            problem = struct('A',[0 1],'b',0,'c',[1;1],'osense',-1,'F',speye(2),'lb',[0;0],'ub',[0;0],'csense','E','vartype',['C';'I'],'x0',[0;0]);
+            problem = struct('A',[0 1],'b',0,'c',[1;1],'osense',-1,'F',speye(2),'lb',[0;0],'ub',[0;0],'csense','E','vartype',['C';'I'],'x0',[0;0],'names',[]);
         end
         try
             % Skip the CLP solver until further developments
@@ -542,9 +542,9 @@
                 %This is the code that actually tests if a solver is working
                 if validationLevel>1
                     %display progress
-                    eval(['solveCobra' solverType '(problem,''printLevel'',3);']);
+                    eval(['solveCobra' problemType '(problem,''printLevel'',3);']);
                 else
-                    eval(['solveCobra' solverType '(problem,''printLevel'',0);']);
+                    eval(['solveCobra' problemType '(problem,''printLevel'',0);']);
                 end
             end
         catch ME
@@ -556,12 +556,12 @@
                 rethrow(ME)
             end
             solverOK = false;
-            eval(['CBT_', solverType, '_SOLVER = oldval;']);
+            eval(['CBT_', problemType, '_SOLVER = oldval;']);
         end
         warning(cwarn)
     else
         % if unvalidated, simply set the solver without testing.
-        eval(['CBT_', solverType, '_SOLVER = solverName;']);
+        eval(['CBT_', problemType, '_SOLVER = solverName;']);
     end
 else
     switch solverName
diff --git a/src/base/solvers/gurobi/setGurobiParam.m b/src/base/solvers/gurobi/setGurobiParam.m
new file mode 100644
index 0000000000..edceb9cf66
--- /dev/null
+++ b/src/base/solvers/gurobi/setGurobiParam.m
@@ -0,0 +1,356 @@
+function gurobiParam = setGurobiParam(param)
+%  The params struct contains Gurobi parameters. A full list may be
+%  found on the Parameter page of the reference manual:
+%  https://www.gurobi.com/documentation/current/refman/parameter_descriptions.html
+%  Parameters must be in TimeLimit not timelimit or timeLimit format, see
+%  below for full list that are eligible to be passed to the solver.
+
+if isfield(param,'timelimit')
+    param.TimeLimit=param.timelimit;
+end
+
+if isfield(param,'multiscale') && param.multiscale==1
+    param.ScaleFlag=0;
+end
+
+if isfield(param,'lifted') && param.lifted==1
+    param.Aggregate=1;
+    % Presolve
+    % Controls the presolve level
+    % Type:	int
+    % Default value:	-1
+    % Minimum value:	-1
+    % Maximum value:	2
+    % Controls the presolve level. A value of -1 corresponds to an automatic setting.
+    % Other options are off (0), conservative (1), or aggressive (2). More aggressive application
+    % of presolve takes more time, but can sometimes lead to a significantly tighter model.
+    param.Presolve=0;
+end
+
+%backward compatibility
+if isfield(param,'method')
+    if isempty(param.method)
+        param = rmfield(param,'method');
+    else
+        if ~isfield(param,[lower(param.problemType) 'method'])
+            param.([lower(param.problemType) 'method'])=param.method;
+        end
+    end
+end
+
+% https://www.gurobi.com/documentation/current/refman/method.html
+% params.method gives the method used to solve continuous models
+% -1=automatic,
+%  0=primal simplex,
+%  1=dual simplex,
+%  2=barrier,
+%  3=concurrent,
+%  4=deterministic concurrent
+% i.e. params.method     = 1;          % use dual simplex method
+if isfield(param,'lpmethod')
+    %gurobiAlgorithms = {'AUTOMATIC','PRIMAL','DUAL','BARRIER','CONCURRENT','CONCURRENT_DETERMINISTIC'};
+    % -1=automatic,
+    % 0=primal simplex,
+    % 1=dual simplex,
+    % 2=barrier,
+    % 3=concurrent,
+    % 4=deterministic concurrent
+    switch param.lpmethod
+        case 'AUTOMATIC'
+            param.Method = -1;
+        case 'PRIMAL'
+            param.Method = 0;
+        case 'DUAL'
+            param.Method = 1;
+        case 'BARRIER'
+            param.Method = 2;
+        case 'CONCURRENT'
+            param.Method = 3;
+        case 'DETERMINISTIC_CONCURRENT'
+            param.Method = 4;
+        otherwise
+            %https://www.gurobi.com/documentation/current/refman/method.html
+            %Concurrent methods aren't available for QP and QCP.
+            warning([param.lpmethod ' is an unrecognised param.qpmethod for gurobi'])
+    end
+    param = rmfield(param,'lpmethod');
+end
+
+% param.qpmethod gives the qpmethod used to solve continuous models
+% -1=automatic,
+%  0=primal simplex,
+%  1=dual simplex,
+%  2=barrier,
+%  3=concurrent,
+%  4=deterministic concurrent
+% i.e. param.qpmethod     = 1;          % use dual simplex method
+if isfield(param,'qpmethod')
+    %gurobi QP algorithms = {'AUTOMATIC','PRIMAL','DUAL','BARRIER'};
+    % -1=automatic,
+    % 0=primal simplex,
+    % 1=dual simplex,
+    % 2=barrier,
+    switch param.qpmethod
+        case 'AUTOMATIC'
+            param.Method = -1;
+        case 'PRIMAL'
+            param.Method = 0;
+        case 'DUAL'
+            param.Method = 1;
+        case 'BARRIER'
+            param.Method = 2;
+        otherwise
+            %https://www.gurobi.com/documentation/current/refman/method.html
+            %Concurrent methods aren't available for QP and QCP.
+            warning([param.qpmethod ' is an unrecognised param.qpmethod for gurobi'])
+    end
+    param = rmfield(param,'qpmethod');
+end
+
+if ~isfield(param,'OutputFlag')
+    switch param.printLevel
+        case 0
+            param.OutputFlag = 0;
+        case 1
+            param.OutputFlag = 0;
+        otherwise
+            % silent
+            param.OutputFlag = 0;
+    end
+end
+
+if ~isfield(param,'DisplayInterval')
+    switch param.printLevel
+        case 0
+            param.DisplayInterval = 1;
+        case 1
+            param.DisplayInterval = 1;
+        otherwise
+            % silent
+            param.DisplayInterval = 1;
+    end
+end
+
+if ~isfield(param,'FeasibilityTol')
+    % Primal feasibility tolerance
+    % Type:	double
+    % Default value:	1e-6
+    % Minimum value:	1e-9
+    % Maximum value:	1e-2
+    % All constraints must be satisfied to a tolerance of FeasibilityTol.
+    % Tightening this tolerance can produce smaller constraint violations, but for
+    % numerically challenging models it can sometimes lead to much larger iteration counts.
+     param.FeasibilityTol = param.feasTol;
+end
+if ~isfield(param,'OptimalityTol')
+    % Dual feasibility tolerance
+    % Type:	double
+    % Default value:	1e-6
+    % Minimum value:	1e-9
+    % Maximum value:	1e-2
+    % Reduced costs must all be smaller than OptimalityTol in the improving direction in order for a model to be declared optimal.
+     param.OptimalityTol = param.optTol;
+end
+
+
+% Permitted parameter fields
+permittedFields = {...
+'AggFill'...
+'Aggregate'...
+'BarConvTol'...
+'BarCorrectors'...
+'BarHomogeneous'...
+'BarIterLimit'...
+'BarOrder'...
+'BarQCPConvTol'...
+'BestBdStop'...
+'BestObjStop'...
+'BQPCuts'...
+'BranchDir'...
+'CliqueCuts'...
+'CloudAccessID'...
+'CloudHost'...
+'CloudSecretKey'...
+'CloudPool'...
+'ComputeServer'...
+'ConcurrentJobs'...
+'ConcurrentMethod'...
+'ConcurrentMIP'...
+'ConcurrentSettings'...
+'CoverCuts'...
+'Crossover'...
+'CrossoverBasis'...
+'CSAPIAccessID'...
+'CSAPISecret'...
+'CSAppName'...
+'CSAuthToken'...
+'CSBatchMode'...
+'CSClientLog'...
+'CSGroup'...
+'CSIdleTimeout'...
+'CSManager'...
+'CSPriority'...
+'CSQueueTimeout'...
+'CSRouter'...
+'CSTLSInsecure'...
+'CutAggPasses'...
+'Cutoff'...
+'CutPasses'...
+'Cuts'...
+'DegenMoves'...
+'Disconnected'...
+'DisplayInterval'...
+'DistributedMIPJobs'...
+'DualReductions'...
+'FeasibilityTol'...
+'FeasRelaxBigM'...
+'FlowCoverCuts'...
+'FlowPathCuts'...
+'FuncPieceError'...
+'FuncPieceLength'...
+'FuncPieceRatio'...
+'FuncPieces'...
+'FuncMaxVal'...
+'FuncNonlinear'...
+'GomoryPasses'...
+'GUBCoverCuts'...
+'Heuristics'...
+'IgnoreNames'...
+'IISMethod'...
+'ImpliedCuts'...
+'ImproveStartGap'...
+'ImproveStartNodes'...
+'ImproveStartTime'...
+'InfProofCuts'...
+'InfUnbdInfo'...
+'InputFile'...
+'IntegralityFocus'...
+'IntFeasTol'...
+'IterationLimit'...
+'JobID'...
+'JSONSolDetail'...
+'LazyConstraints'...
+'LicenseID'...
+'LiftProjectCuts'...
+'LPWarmStart'...
+'LogFile'...
+'LogToConsole'...
+'MarkowitzTol'...
+'MemLimit'...
+'Method'...
+'MinRelNodes'...
+'MIPFocus'...
+'MIPGap'...
+'MIPGapAbs'...
+'MIPSepCuts'...
+'MIQCPMethod'...
+'MIRCuts'...
+'MixingCuts'...
+'ModKCuts'...
+'MultiObjMethod'...
+'MultiObjPre'...
+'MultiObjSettings'...
+'NetworkAlg'...
+'NetworkCuts'...
+'NLPHeur'...
+'NodefileDir'...
+'NodefileStart'...
+'NodeLimit'...
+'NodeMethod'...
+'NonConvex'...
+'NoRelHeurTime'...
+'NoRelHeurWork'...
+'NormAdjust'...
+'NumericFocus'...
+'OBBT'...
+'ObjNumber'...
+'ObjScale'...
+'OptimalityTol'...
+'OutputFlag'...
+'PartitionPlace'...
+'PerturbValue'...
+'PoolGap'...
+'PoolGapAbs'...
+'PoolSearchMode'...
+'PoolSolutions'...
+'PreCrush'...
+'PreDepRow'...
+'PreDual'...
+'PreMIQCPForm'...
+'PrePasses'...
+'PreQLinearize'...
+'Presolve'...
+'PreSOS1BigM'...
+'PreSOS1Encoding'...
+'PreSOS2BigM'...
+'PreSOS2Encoding'...
+'PreSparsify'...
+'ProjImpliedCuts'...
+'PSDCuts'...
+'PSDTol'...
+'PumpPasses'...
+'QCPDual'...
+'Quad'...
+'Record'...
+'ResultFile'...
+'RINS'...
+'RelaxLiftCuts'...
+'RLTCuts'...
+'ScaleFlag'...
+'ScenarioNumber'...
+'Seed'...
+'ServerPassword'...
+'ServerTimeout'...
+'Sifting'...
+'SiftMethod'...
+'SimplexPricing'...
+'SoftMemLimit'...
+'SolutionLimit'...
+'SolutionTarget'...
+'SolFiles'...
+'SolutionNumber'...
+'StartNodeLimit'...
+'StartNumber'...
+'StrongCGCuts'...
+'SubMIPCuts'...
+'SubMIPNodes'...
+'Symmetry'...
+'Threads'...
+'TimeLimit'...
+'TokenServer'...
+'TSPort'...
+'TuneBaseSettings'...
+'TuneCleanup'...
+'TuneCriterion'...
+'TuneDynamicJobs'...
+'TuneJobs'...
+'TuneMetric'...
+'TuneOutput'...
+'TuneResults'...
+'TuneTargetMIPGap'...
+'TuneTargetTime'...
+'TuneTimeLimit'...
+'TuneTrials'...
+'TuneUseFilename'...
+'UpdateMode'...
+'UserName'...
+'VarBranch'...
+'WLSAccessID'...
+'WLSSecret'...
+'WLSToken'...
+'WLSTokenDuration'...
+'WLSTokenRefresh'...
+'WorkerPassword'...
+'WorkerPool'...
+'WorkLimit'...
+'ZeroHalfCuts'...
+'ZeroObjNodes'};
+
+fields = fieldnames(param);
+fieldsToRemove = setdiff(fields, permittedFields);
+%remove fields that should not be present
+gurobiParam = rmfield(param, fieldsToRemove);
+
+if isfield(gurobiParam,'logFile') && isempty(gurobiParam.logFile)
+    gurobiParam = rmfield(gurobiParam,'logFile');
+end
diff --git a/src/base/solvers/init/configEnvVars.m b/src/base/solvers/init/configEnvVars.m
index 89c40b6c92..5d9eb639a2 100644
--- a/src/base/solvers/init/configEnvVars.m
+++ b/src/base/solvers/init/configEnvVars.m
@@ -191,7 +191,7 @@
                 end
 
                 % Method 1: check if the solver is already on the MATLAB path
-                isOnPath = ~isempty(strfind(lower(path), lower(possibleDir)));
+                isOnPath = contains(lower(path), lower(possibleDir));
 
                 % find the index of the most recently added solver path
                 tmp = path;
diff --git a/src/base/solvers/entropicFBA/mosekParamStrip.m b/src/base/solvers/mosek/mosekParamStrip.m
similarity index 100%
rename from src/base/solvers/entropicFBA/mosekParamStrip.m
rename to src/base/solvers/mosek/mosekParamStrip.m
diff --git a/src/base/solvers/msk/parseMskResult.m b/src/base/solvers/mosek/parseMskResult.m
similarity index 92%
rename from src/base/solvers/msk/parseMskResult.m
rename to src/base/solvers/mosek/parseMskResult.m
index 0a61e1081b..05eaf0d088 100644
--- a/src/base/solvers/msk/parseMskResult.m
+++ b/src/base/solvers/mosek/parseMskResult.m
@@ -1,4 +1,4 @@
-function [stat,origStat,x,y,yl,yu,z,zl,zu,k,basis,pobjval,dobjval] = parseMskResult(res,solverOnlyParams,printLevel)
+function [stat,origStat,x,y,yl,yu,z,zl,zu,s,basis,pobjval,dobjval] = parseMskResult(res)
 %parse the res structure returned from mosek
 % INPUTS:
 %  res:        mosek results structure returned by mosekopt
@@ -49,17 +49,11 @@
 z = [];
 zl = [];
 zu = [];
-k = [];
+s = [];
 basis = [];
 pobjval =[];
 dobjval =[];
 
-if ~exist('printLevel','var')
-    printLevel = 0;
-end
-if ~exist('solverOnlyParams','var')
-    solverOnlyParams = struct();
-end
 
 % prosta (string) – Problem status (prosta).
 % prosta
@@ -126,7 +120,7 @@
 % The primal solution is integer optimal.
 
 % https://docs.mosek.com/latest/toolbox/accessing-solution.html
-accessSolution=[];
+accessSolution='';
 if isfield(res, 'sol')
     if isfield(res.sol,'itr') && isfield(res.sol,'bas')
         if  any(strcmp(res.sol.bas.solsta,{'OPTIMAL','MSK_SOL_STA_OPTIMAL','MSK_SOL_STA_NEAR_OPTIMAL'})) && any(strcmp(res.sol.itr.solsta,{'UNKNOWN'}))
@@ -178,7 +172,7 @@
                 zu=res.sol.itr.sux; %dual to   x <= bux
                 if isfield(res.sol.itr,'doty')
                     % Dual variables to affine conic constraints
-                    k = res.sol.itr.doty;
+                    s = res.sol.itr.doty;
                 end
                 pobjval = res.sol.itr.pobjval;
                 dobjval = res.sol.itr.dobjval;
@@ -204,19 +198,22 @@
                 zu=res.sol.bas.sux; %dual to   x <= bux
                 if isfield(res.sol.bas,'s')
                     % Dual variables to affine conic constraints
-                    k = res.sol.bas.s;
+                    s = res.sol.bas.s;
                 end
 
                 %https://docs.mosek.com/10.0/toolbox/advanced-hotstart.html
-                bas.skc = res.sol.bas.skc;
-                bas.skx = res.sol.bas.skx;
-                bas.xc = res.sol.bas.xc;
-                bas.xx = res.sol.bas.xx;
+                basis.skc = res.sol.bas.skc;
+                basis.skx = res.sol.bas.skx;
+                basis.xc = res.sol.bas.xc;
+                basis.xx = res.sol.bas.xx;
                 pobjval = res.sol.bas.pobjval;
                 dobjval = res.sol.bas.dobjval;
             otherwise
                 accessSolution = 'dontAccess';
         end
+    otherwise
+        accessSolution = 'dontAccess';
+        origStat = -1;
 end
 
 if strcmp(accessSolution,'dontAccess')
@@ -227,11 +224,11 @@
         case {'DUAL_INFEASIBLE_CER','MSK_SOL_STA_DUAL_INFEAS_CER','MSK_SOL_STA_NEAR_DUAL_INFEAS_CER'}
             stat=2; % Unbounded solution
             origStat = [origStat ' & ' res.rcodestr];
-        case {'UNKNOWN','PRIM_ILLPOSED_CER','DUAL_ILLPOSED_CER','PRIM_FEAS','DUAL_FEAS','PRIM_AND_DUAL_FEAS'}
+        case {'UNKNOWN','PRIM_ILLPOSED_CER','DUAL_ILLPOSED_CER','PRIM_FEAS','DUAL_FEAS','PRIM_AND_DUAL_FEAS','DUAL_FEASIBLE'}
             stat=-1; %some other problem
             origStat = [origStat ' & ' res.rcodestr];
         otherwise
-            warning(['Unrecognised res.sol.bas.solsta: ' origStat])
+            warning(['Unrecognised res.sol.bas.solsta or res.sol.itr.solsta: ' origStat])
             stat=-1; %some other problem
             fprintf('%s\n',res.rcode)
             fprintf('%s\n',res.rmsg)
@@ -241,3 +238,6 @@
             end
     end
 end
+
+
+% https://themosekblog.blogspot.com/2014/06/what-if-solver-stall.html
\ No newline at end of file
diff --git a/src/base/solvers/mosek/setMosekParam.m b/src/base/solvers/mosek/setMosekParam.m
new file mode 100644
index 0000000000..8b786f0b11
--- /dev/null
+++ b/src/base/solvers/mosek/setMosekParam.m
@@ -0,0 +1,384 @@
+function [cmd,mosekParam] = setMosekParam(param)
+% set mosek parameters from param fields
+% strip any non mosek compatible fields from param and return it as
+% mosekParam
+
+
+%tests if solver correctly interfaced and licence running
+if param.printLevel>1 || param.debug
+    [~, res] = mosekopt('symbcon');
+else
+    [~, res] = mosekopt('symbcon echo(0)');
+end
+
+% only set the print level if not already set via param structure
+if ~isfield(param, 'MSK_IPAR_LOG')
+    % Controls the amount of log information.
+    % The value 0 implies that all log information is suppressed.
+    % A higher level implies that more information is logged.
+    switch param.printLevel
+        case 0
+            echolev = 0;
+        case 1
+            echolev = 3;
+        case 2
+            param.MSK_IPAR_WRITE_DATA_PARAM='MSK_ON';
+            param.MSK_IPAR_LOG_INTPNT = 1;
+            param.MSK_IPAR_LOG_SIM = 1;
+            %MSK_IPAR_LOG_PRESOLVE
+            % Description:Controls amount of output printed by the presolve procedure. A higher level implies that more information is logged.
+            % Possible Values:Any number between 0 and +inf.
+            % Default value:1
+            param.MSK_IPAR_LOG_PRESOLVE=10;
+
+            %MSK_IPAR_LOG_INTPNT
+            % Controls amount of output printed printed by the interior-point optimizer.
+            %A higher level implies that more information is logged.
+            % Possible Values: Any number between 0 and +inf.
+            % Default value: 4
+            if ~isfield(param,'MSK_IPAR_LOG_INTPNT')
+                param.MSK_IPAR_LOG_INTPNT=5;
+            end
+
+            %infesibility report
+            % MSK_IPAR_INFEAS_REPORT_AUTO
+            %Controls the amount of information presented in an infeasibility report.
+            % Possible Values:
+            %     MSK_ON
+            %         Switch the option on.
+            %     MSK_OFF
+            %         Switch the option off.
+            % Default value:
+            %     MSK_OFF
+            if ~isfield(param,'MSK_IPAR_INFEAS_REPORT_AUTO')
+                param.MSK_IPAR_INFEAS_REPORT_AUTO='MSK_ON';
+            end
+
+            % MSK_IPAR_INFEAS_REPORT_LEVEL
+            % Controls the amount of information presented in an infeasibility report. Higher values imply more information.
+            % Possible Values:Any number between 0 and +inf.
+            % Default value: 1
+            %Higher values imply more information.
+            if ~isfield(param,'MSK_IPAR_INFEAS_REPORT_LEVEL')
+                param.MSK_IPAR_INFEAS_REPORT_LEVEL=1;
+            end
+
+            echolev = 3;
+        otherwise
+            echolev = 0;
+    end
+    if echolev == 0 && ~param.debug
+        param.MSK_IPAR_LOG = 0;
+        cmd = ['minimize echo(' int2str(echolev) ')'];
+    else
+        cmd = 'minimize';
+    end
+end
+
+
+if ~isfield(param, 'MSK_DPAR_OPTIMIZER_MAX_TIME') && isfield(param,'timelimit')
+    % MSK_DPAR_OPTIMIZER_MAX_TIME
+    % Maximum amount of time the optimizer is allowed to spent on the optimization. A negative number means infinity.
+    % Default
+    % -1.0
+    % Accepted
+    % [-inf; +inf]
+    % Example
+    % param.MSK_DPAR_OPTIMIZER_MAX_TIME = -1.0
+    % Groups
+    % Termination criteria
+    param.MSK_DPAR_OPTIMIZER_MAX_TIME = param.timelimit;
+end
+
+
+
+if ~isfield(param, 'MSK_DPAR_INTPNT_TOL_PFEAS')
+    % Primal feasibility tolerance used by the interior-point optimizer for linear problems.
+    % Default
+    % 1.0e-8
+    % Accepted
+    % [0.0; 1.0]
+    % Example
+    % param.MSK_DPAR_INTPNT_TOL_PFEAS = 1.0e-8
+    % Groups
+    % Interior-point method, Termination criteria
+    param.MSK_DPAR_INTPNT_TOL_PFEAS=param.feasTol;
+end
+
+
+if ~isfield(param,'MSK_DPAR_INTPNT_QO_TOL_PFEAS')
+    % Primal feasibility tolerance used by the interior-point optimizer for quadratic problems.
+    % Default
+    % 1.0e-8
+    % Accepted
+    % [0.0; 1.0]
+    % Example
+    % param.MSK_DPAR_INTPNT_QO_TOL_PFEAS = 1.0e-8
+    % See also
+    % MSK_DPAR_INTPNT_QO_TOL_NEAR_REL
+    % Groups
+    % Interior-point method, Termination criteria
+    param.MSK_DPAR_INTPNT_QO_TOL_PFEAS=param.feasTol;
+end
+
+if ~isfield(param, 'MSK_DPAR_INTPNT_CO_TOL_PFEAS')
+    % Primal feasibility tolerance used by the interior-point optimizer for conic problems.
+    % Default
+    % 1.0e-8
+    % Accepted
+    % [0.0; 1.0]
+    % Example
+    % param.MSK_DPAR_INTPNT_CO_TOL_PFEAS = 1.0e-8
+    % See also
+    % MSK_DPAR_INTPNT_CO_TOL_NEAR_REL
+    % Groups
+    % Interior-point method, Termination criteria, Conic interior-point method
+    param.MSK_DPAR_INTPNT_CO_TOL_PFEAS=param.feasTol;
+end
+
+if ~isfield(param, 'MSK_DPAR_INTPNT_TOL_DFEAS')
+    % MSK_DPAR_INTPNT_TOL_DFEAS
+    % Dual feasibility tolerance used by the interior-point optimizer for linear problems.
+    % Default
+    % 1.0e-8
+    % Accepted
+    % [0.0; 1.0]
+    % Example
+    % param.MSK_DPAR_INTPNT_TOL_DFEAS = 1.0e-8
+    % Groups
+    % Interior-point method, Termination criteria
+    param.MSK_DPAR_INTPNT_TOL_DFEAS=param.optTol;
+end
+
+if ~isfield(param, 'MSK_DPAR_INTPNT_QO_TOL_DFEAS')
+    % Dual feasibility tolerance used by the interior-point optimizer for quadratic problems.
+    % Default
+    % 1.0e-8
+    % Accepted
+    % [0.0; 1.0]
+    % Example
+    % param.MSK_DPAR_INTPNT_QO_TOL_DFEAS = 1.0e-8
+    % See also
+    % MSK_DPAR_INTPNT_QO_TOL_NEAR_REL
+    % Groups
+    % Interior-point method, Termination criteria
+    param.MSK_DPAR_INTPNT_QO_TOL_DFEAS=param.optTol;
+end
+
+if ~isfield(param, 'MSK_DPAR_INTPNT_CO_TOL_DFEAS')
+    % Dual feasibility tolerance used by the interior-point optimizer for linear problems.
+    % Default
+    % 1.0e-8
+    % Accepted
+    % [0.0; 1.0]
+    % Example
+    % param.MSK_DPAR_INTPNT_TOL_DFEAS = 1.0e-8
+    % Groups
+    % Interior-point method, Termination criteria
+    param.MSK_DPAR_INTPNT_CO_TOL_DFEAS=param.optTol;
+end
+
+if isfield(param,'lifted') && param.lifted==1
+    % Controls the maximum amount of fill-in that can be created by one pivot in the elimination phase of the presolve.
+    % A negative value means the parameter value is selected automatically.
+    % Default-1
+    % Accepted [-inf; +inf]
+    % Example param.MSK_IPAR_PRESOLVE_ELIMINATOR_MAX_FILL = -1
+    if ~isfield(param,'MSK_IPAR_PRESOLVE_ELIMINATOR_MAX_NUM_TRIES')
+        param.MSK_IPAR_PRESOLVE_ELIMINATOR_MAX_NUM_TRIES = 0;
+    end
+end
+
+%turn on multiscale if infeasibilities after unscaling
+if isfield(param,'multiscale') && param.multiscale==1 && param.lifted==0
+    % Controls whether whether a new experimental linear dependency checker is employed.
+    % Default
+    % "OFF"
+    % Accepted
+    % "ON", "OFF"
+    % Example
+    % param.MSK_IPAR_PRESOLVE_LINDEP_NEW = 'MSK_OFF'
+    if ~isfield(param,'MSK_IPAR_PRESOLVE_LINDEP_NEW')
+        param.MSK_IPAR_PRESOLVE_LINDEP_NEW = 'MSK_OFF';
+    end
+    
+    % MSK_IPAR_INTPNT_SCALING
+    % Controls how the problem is scaled before the interior-point optimizer is used.
+    % Default
+    % "FREE"
+    % Accepted
+    % "FREE", "NONE"
+    % param..MSK_IPAR_INTPNT_SCALING = 'MSK_SCALING_FREE';
+    if ~isfield(param,'MSK_IPAR_INTPNT_SCALING')
+        param.MSK_IPAR_INTPNT_SCALING='MSK_SCALING_NONE';
+    end
+    % MSK_IPAR_SIM_SCALING
+    % Controls how much effort is used in scaling the problem before a simplex optimizer is used.
+    % Default
+    % "FREE"
+    % Accepted
+    % "FREE", "NONE"
+    % Example
+    % param.MSK_IPAR_SIM_SCALING = 'MSK_SCALING_FREE'
+    if ~isfield(param,'MSK_IPAR_SIM_SCALING')
+        param.MSK_IPAR_SIM_SCALING='MSK_SCALING_NONE';
+    end
+end
+
+if isfield(param,'debug') && param.debug==1
+    % https://docs.mosek.com/latest/rmosek/debugging-infeas.html
+    % Controls whether an infeasibility report is automatically produced after the optimization if the problem is primal or dual infeasible.
+    param.MSK_IPAR_INFEAS_REPORT_AUTO='MSK_ON';
+end
+
+if isfield(param,'strict')
+    % MSK_IPAR_BI_IGNORE_MAX_ITER
+    % If the parameter MSK_IPAR_INTPNT_BASIS has the value MSK_BI_NO_ERROR and
+    % the interior-point optimizer has terminated due to maximum number of
+    % iterations, then basis identification is performed if this parameter has
+    % the value MSK_ON.
+    % Possible Values:
+    %     MSK_ON        Switch the option on.
+    %     MSK_OFF       Switch the option off.
+    % Default value:
+    %     MSK_OFF
+    if ~isfield(param,'MSK_IPAR_BI_IGNORE_MAX_ITER')
+        param.MSK_IPAR_BI_IGNORE_MAX_ITER='MSK_OFF';
+    end
+
+    %%%%%%%%%%% Solution Approach
+    % MSK_IPAR_INTPNT_SOLVE_FORM
+    % Controls whether the primal or the dual problem is solved.
+    % Possible Values:
+    %     MSK_SOLVE_PRIMAL
+    %         The optimizer should solve the primal problem.
+    %     MSK_SOLVE_DUAL
+    %         The optimizer should solve the dual problem.
+    %     MSK_SOLVE_FREE
+    %         The optimizer is free to solve either the primal or the dual problem.
+    % Default value:MSK_SOLVE_FREE
+    if ~isfield(param,'MSK_IPAR_INTPNT_SOLVE_FORM')
+        param.MSK_IPAR_INTPNT_SOLVE_FORM='MSK_SOLVE_FREE';
+        %param.MSK_IPAR_INTPNT_SOLVE_FORM='MSK_SOLVE_PRIMAL';
+    end
+
+    %%%%%%% Infeasibility
+    % MSK_DPAR_INTPNT_TOL_INFEAS
+    % Controls when the optimizer declares the model primal or dual infeasible.
+    % A small number means the optimizer gets more conservative about declaring the model infeasible.
+    % Possible Values:Any number between 0.0 and 1.0.
+    % Default value: 1.0e-8
+    if ~isfield(param,'MSK_DPAR_INTPNT_TOL_INFEAS')
+        % param.MSK_DPAR_INTPNT_TOL_INFEAS=1e-10;
+        param.MSK_DPAR_INTPNT_TOL_INFEAS=1e-8;
+    end
+end
+
+%backward compatibility
+if isfield(param,'method')
+    if isempty(param.method)
+        param = rmfield(param,'method');
+    else
+        if ~isfield(param,[lower(param.problemType) 'method'])
+            param.([lower(param.problemType) 'method'])=param.method;
+        end
+    end
+end
+
+switch param.problemType
+    case {'LP'}
+        if isfield(param,'lpmethod')
+            if contains(param.lpmethod,'MSK_OPTIMIZER_')
+                param.MSK_IPAR_OPTIMIZER=param.lpmethod;
+            else
+                param.MSK_IPAR_OPTIMIZER=['MSK_OPTIMIZER_' param.lpmethod];
+            end
+        end
+    case {'QP'}
+        if isfield(param,'qpmethod')
+            if contains(param.qpmethod,'MSK_OPTIMIZER_')
+                param.MSK_IPAR_OPTIMIZER=param.qpmethod;
+            else
+                param.MSK_IPAR_OPTIMIZER=['MSK_OPTIMIZER_' param.qpmethod];
+            end
+        end
+    case {'CLP'}
+        if isfield(param,'clpmethod')
+            if contains(param.qpmethod,'MSK_OPTIMIZER_')
+                param.MSK_IPAR_OPTIMIZER=param.clpmethod;
+            else
+                param.MSK_IPAR_OPTIMIZER=['MSK_OPTIMIZER_' param.clpmethod];
+            end
+        end
+
+        % MSK_IPAR_INTPNT_REGULARIZATION_USE
+        % Description:Controls whether regularization is allowed.
+        % Possible Values: MSK_ON    Switch the option on.
+        %                  MSK_OFF   Switch the option off.
+        % Default value:   MSK_ON
+        if ~isfield(param,'MSK_IPAR_INTPNT_REGULARIZATION_USE')
+            param.MSK_IPAR_INTPNT_REGULARIZATION_USE='MSK_ON';
+        end
+
+    case {'EP'}
+        if isfield(param,'epmethod')
+            if contains(param.epmethod,'MSK_OPTIMIZER_')
+                param.MSK_IPAR_OPTIMIZER=param.epmethod;
+            else
+                param.MSK_IPAR_OPTIMIZER=['MSK_OPTIMIZER_' param.epmethod];
+            end
+        end
+        % MSK_IPAR_INTPNT_REGULARIZATION_USE
+        % Description:Controls whether regularization is allowed.
+        % Possible Values: MSK_ON    Switch the option on.
+        %                  MSK_OFF   Switch the option off.
+        % Default value:   MSK_ON
+        if ~isfield(param,'MSK_IPAR_INTPNT_REGULARIZATION_USE')
+            param.MSK_IPAR_INTPNT_REGULARIZATION_USE='MSK_ON';
+        end
+
+
+
+        % MSK_IPAR_INTPNT_MAX_ITERATIONS
+        % Controls the maximum number of iterations allowed in the interior-point optimizer.
+        % Possible Values:Any number between 0 and +inf.
+        % Default value: 400
+        if ~isfield(param,'MSK_IPAR_INTPNT_MAX_ITERATIONS')
+            param.MSK_IPAR_INTPNT_MAX_ITERATIONS=400;
+        end
+
+        %%%%%%%%%%%%% NONLINEAR SOLVER INTEGER PARAM %%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+    case {'VK'}
+
+        % MSK_IPAR_SIM_SCALING_METHOD
+        % Controls how the problem is scaled before a simplex optimizer is used.
+        % Default
+        % "POW2"
+        % Accepted
+        % "POW2", "FREE"
+        % Example
+        % param.MSK_IPAR_SIM_SCALING_METHOD = 'MSK_SCALING_METHOD_POW2'
+        % param.MSK_IPAR_SIM_SCALING_METHOD='MSK_SCALING_METHOD_FREE';
+end
+
+
+if ~isfield(param,'MSK_IPAR_LOG_FEAS_REPAIR') && isfield(param,'repairInfeasibility')
+    % MSK_IPAR_LOG_FEAS_REPAIR
+    % Controls the amount of output printed when performing feasibility repair. A value higher than one means extensive logging.
+    % Default
+    % 1
+    % Accepted
+    % [0; +inf]
+    % Example
+    % MSK_putintparam(task, MSK_IPAR_LOG_FEAS_REPAIR, 1)
+    % Groups
+    % Output information, Logging
+    param.MSK_IPAR_LOG_FEAS_REPAIR = param.repairInfeasibility;
+end
+
+% Remove outer function specific parameters to avoid crashing solver interfaces
+mosekParam = mosekParamStrip(param);
+if 0
+    disp(mosekParam)
+end
diff --git a/src/base/solvers/param/getCobraSolverParams.m b/src/base/solvers/param/getCobraSolverParams.m
index 96f94931bf..92daa82e30 100644
--- a/src/base/solvers/param/getCobraSolverParams.m
+++ b/src/base/solvers/param/getCobraSolverParams.m
@@ -1,6 +1,6 @@
 function varargout = getCobraSolverParams(problemType, paramNames, paramStructure)
 % This function gets the specified paramStructure in `paramNames` from
-% paramStructure, the global cobra paramters variable or default values set within
+% paramStructure, the global cobra parameters variable or default values set within
 % this script. 
 %
 % It will use values with the following priority
@@ -43,21 +43,47 @@
 %       - Richard Que (12/01/2009)
 %       - Ronan (16/07/2013) default MPS paramStructure are no longer global variables
 
-if nargin < 2
-    error('getCobraSolverParams: No paramStructure specified')
-end
-
-if nargin < 3
-    paramStructure = [];
+if ~exist('paramNames','var') || isempty(paramNames)
+    %get the names of all the parameters
+    defaultParamNames = getCobraSolverParamsOptionsForType(problemType);
+else
+    if ~iscell(paramNames)
+        paramNames = {paramNames};
+    end
 end
-if ~isempty(paramStructure)
-    paramStructureOut = paramStructure;
+if exist('paramStructure','var')
+    overRideDefaults = 1;
+    %overide defaults unless hack is asking for defaults
+    if isstruct(paramStructure)
+        paramStructureOut = paramStructure;
+    else
+        if ischar(paramStructure) && strcmp(paramStructure,'default')
+           overRideDefaults=1;
+        end
+    end
+else
+    %the last argument will be the structure with defaults
+    overRideDefaults=0;
+    paramStructure=[];
 end
+
 % Persistence will make names specific to one type of solver.
 % Default Values
 % For descriptions of the different settings please have a look at 
 % getCobraSolverParamsOptionsForType
-valDef.minNorm = 0;
+
+% %These default tolerances are based on the default values for the Gurobi LP
+% %solver. Do not change them without first consulting with other developers.
+% %https://www.gurobi.com/documentation/9.0/refman/parameters.html
+% % (primal) feasibility tolerance
+% changeCobraSolverParams('LP', 'feasTol', 1e-6);
+% % (dual) optimality tolerance
+% changeCobraSolverParams('LP', 'optTol', 1e-6);
+% 
+% % (primal) feasibility tolerance
+% changeCobraSolverParams('EP', 'feasTol', 1e-8);
+% % (dual) optimality tolerance
+% changeCobraSolverParams('EP', 'optTol', 1e-12);
 
 %These default tolerances are based on the default values for the Gurobi LP solver
 %https://www.gurobi.com/documentation/9.0/refman/paramStructure.html
@@ -76,20 +102,12 @@
 valDef.timeLimit = 1e36; 
 valDef.iterationLimit = 1000;
 
-%valDef.logFile = ['Cobra' problemType 'Solver.log'];
 valDef.logFile = []; %log file should be empty to avoid creating it by default
 valDef.saveInput = [];
-valDef.PbName = [problemType 'problem'];
+valDef.problemType = problemType;
 valDef.debug = 0;
 valDef.lifting = 0;
-   
-valDef.method = -1; 
-
-% CPLEX paramStructure
-valDef.DATACHECK = 1;
-valDef.DEPIND = 1;
-valDef.checkNaN = 0;
-valDef.warning = 0;
+valDef.multiscale=0;  % true if problem is multiscale
 
 % tolerances
 valDef.intTol = 1e-12;
@@ -97,10 +115,6 @@
 valDef.absMipGapTol = 1e-12;
 valDef.NUMERICALEMPHASIS = 1;    
 
-if ~iscell(paramNames)
-    paramNames = {paramNames};
-end
-
 switch problemType
     case 'LP'
         global CBT_LP_PARAMS
@@ -117,9 +131,13 @@
         valDef.feasTol = 1e-6; % (primal) feasibility tolerance
         valDef.optTol = 1e-6;  % (dual) optimality tolerance
         valDef.solver='mosek';
+        
     case 'CLP'
+        % This is never used elsewhere except for parameter setting loop
+        % for backward compatibility
         global CBT_CLP_PARAMS
         parametersGlobal = CBT_CLP_PARAMS;
+
         valDef.feasTol = 1e-6; % (primal) feasibility tolerance
         valDef.optTol = 1e-6;  % (dual) optimality tolerance
         valDef.solver='mosek';
@@ -134,27 +152,44 @@
         return;
 end
 
+if exist('defaultParamNames','var')
+    % set parameter structure values to default, unless they are overridden by paramStructure
+    for i=1:length(defaultParamNames)
+        if isfield(valDef,defaultParamNames{i})
+            %structure of default parameter values
+            paramStructureOut.(defaultParamNames{i})=valDef.(defaultParamNames{i});
+        end
+    end
+end
+
 varargout = cell(1, numel(paramNames));
 paramNames = columnVector(paramNames);
 for i=1:length(paramNames)
-    % set values to default
+    % first set each value to default
     if isfield(valDef,paramNames{i})
+        %list of default parameter values 
         varargout{i} = valDef.(paramNames{i});
+        %structure of default parameter values
+        paramStructureOut.(paramNames{i})= valDef.(paramNames{i});
     end
-    if ~strcmp(paramStructure,'default') % skip of using default values
-        % set values to global values
+
+    if overRideDefaults
+        % second set values to global values
         if isfield(parametersGlobal,paramNames{i})
-           varargout{i} = parametersGlobal.(paramNames{i});
+            %list of default parameter values 
+            varargout{i} = parametersGlobal.(paramNames{i});
+            %structure of default parameter values
+            paramStructureOut.(paramNames{i})=parametersGlobal.(paramNames{i});
         end
-        % set values to specified values
+        % third set values to specified values in paramStructure
         if isfield(paramStructure,paramNames{i})
-           varargout{i} = paramStructure.(paramNames{i});
-           %remove this field so paramStructureOut can be passed directly to the solver, without extra fields
-           paramStructureOut = rmfield(paramStructureOut,paramNames{i});
+            varargout{i} = paramStructure.(paramNames{i});
+            %structure of default parameter values
+            paramStructureOut.(paramNames{i})=paramStructure.(paramNames{i});
         end
     end
 end
-if ~isempty(paramStructure)
+if overRideDefaults
     varargout{length(paramNames)+1}=paramStructureOut;
 end
 
diff --git a/src/base/solvers/param/getCobraSolverParamsOptionsForType.m b/src/base/solvers/param/getCobraSolverParamsOptionsForType.m
index d228910f21..570eac374c 100644
--- a/src/base/solvers/param/getCobraSolverParamsOptionsForType.m
+++ b/src/base/solvers/param/getCobraSolverParamsOptionsForType.m
@@ -7,24 +7,26 @@
 %
 % INPUT:
 %    problemType :      One of the problem types available in the COBRA
-%                       Toolbox ('LP','QP','MILP','MIQP','NLP')
+%                       Toolbox ('LP','QP','EP','CLP','MILP','MIQP','NLP')
 %
 % OUPTUT:
 %    paramNames:        Cell array of names of parameters that can be set
 %                       for each problem type, independent of the specific
 %                       solver being used.
 
-if iscell(problemType )
+if iscell(problemType)
     paramNames = {};
-    for j = 1:numel(problemType )
-        paramNames = [paramNames, getCobraSolverParamsOptionsForType(problemType {j})];
+    for j = 1:numel(problemType)
+        paramNames = [paramNames, getCobraSolverParamsOptionsForType(problemType{j})];
     end
     paramNames = unique(paramNames);
     return
 end
 switch problemType 
     case 'LP'
-        paramNames = {'verify',...          % verify that it is a suitable  LP problem
+        paramNames = {'multiscale'...       % true if problem is multiscale
+                      'problemType'...      % problem type
+                      'verify',...          % verify that it is a suitable  LP problem
                       'minNorm', ...        % type of normalization used.
                       'printLevel', ...     % print Level
                       'primalOnly', ...     % only solve for primal
@@ -34,11 +36,11 @@
                       'solver', ...         % solver to use (overriding set solver)
                       'debug', ...          % run debgugging code
                       'logFile', ...        % file (location) to write logs to
-                      'lifting', ...        % whether to lift a problem
-                      'method'};            % solver method: -1 = automatic, 0 = primal simplex, 1 = dual simplex, 2 = barrier, 3 = concurrent, 4 = deterministic concurrent, 5 = Network Solver(if supported by the solver)
+                              };            % 
    case 'QP'
-        paramNames = {'verify',...          % verify that it is a suitable  QP problem
-                      'method', ...         % solver method: -1 = automatic, 0 = primal simplex, 1 = dual simplex, 2 = barrier, 3 = concurrent, 4 = deterministic concurrent, 5 = Network Solver(if supported by the solver)
+        paramNames = {'multiscale'...       % true if problem is multiscale
+                      'problemType'...      % problem type
+                      'verify',...          % verify that it is a suitable  QP problem
                       'printLevel', ...     % print level
                       'saveInput', ...      % save the input to a file (specified)
                       'debug', ...          % run debgugging code
@@ -48,15 +50,18 @@
                       'solver'};            % the solver to use
 
    case 'EP'
-        paramNames = {'verify',...          % verify that it is a suitable  EP problem
-                      'method', ...         % solver method: -1 = automatic, 0 = primal simplex, 1 = dual simplex, 2 = barrier, 3 = concurrent, 4 = deterministic concurrent, 5 = Network Solver(if supported by the solver)
+        paramNames = {'multiscale'...       % true if problem is multiscale
+                      'problemType'...      % problem type
+                      'verify',...          % verify that it is a suitable  EP problem
                       'printLevel', ...     % print level
                       'debug', ...          % run debgugging code
                       'feasTol',...         % feasibility tolerance
                       'optTol',...          % optimality tolerance
                       'solver'};            % the solver to use
    case 'CLP'
-        paramNames = {'verify',...          % verify that it is a suitable  CLP problem
+        paramNames = {'multiscale'...       % true if problem is multiscale
+                      'problemType'...      % problem type
+                      'verify',...          % verify that it is a suitable  CLP problem
                       'printLevel', ...     % print level
                       'debug', ...          % run debgugging code
                       'feasTol',...         % feasibility tolerance
@@ -65,7 +70,8 @@
                       };            
                                  
     case 'MILP'
-        paramNames = {'intTol', ...         % integer tolerance (accepted derivation from integer numbers)
+        paramNames = {'problemType'...      % problem type
+                      'intTol', ...         % integer tolerance (accepted derivation from integer numbers)
                       'relMipGapTol', ...   % relative MIP Gap tolerance
                       'absMipGapTol', ...   % absolute MIP Gap tolerance
                       'timeLimit', ...      % maximum time before stopping computation (if supported by the solver)
@@ -78,7 +84,8 @@
                       'debug'};             % run debgugging code
 
     case 'MIQP'
-        paramNames = {'timeLimit', ...      % maximum time before stopping computation (if supported by the solver)
+        paramNames = {'problemType'...      % problem type
+                      'timeLimit', ...      % maximum time before stopping computation (if supported by the solver)
                       'method', ...         % solver method: -1 = automatic, 0 = primal simplex, 1 = dual simplex, 2 = barrier, 3 = concurrent, 4 = deterministic concurrent, 5 = Network Solver(if supported by the solver)
                       'feasTol',...         % feasibility tolerance
                       'optTol',...          % optimality tolerance
@@ -91,7 +98,8 @@
                       'solver'};            % the solver to use
 
     case 'NLP'
-        paramNames = {'warning', ...        % whether to display warnings
+        paramNames = {'problemType'...      % problem type
+                      'warning', ...        % whether to display warnings
                       'checkNaN', ...       % check for NaN solutions
                       'PbName', ...         % name of the problem
                       'iterationLimit', ... % maximum number of iterations before stopping computation (if supported by the solver)
diff --git a/src/base/solvers/param/mergeCobraParams.m b/src/base/solvers/param/mergeCobraParams.m
new file mode 100644
index 0000000000..6e3e991a57
--- /dev/null
+++ b/src/base/solvers/param/mergeCobraParams.m
@@ -0,0 +1,11 @@
+function param = mergeCobraParams(param1,param2)
+% Merge the structures with preference to param1 fields over param2 fields
+
+param = param2;  % Start with the first structure
+fields = fieldnames(param1);  % Get the field names of the second structure
+for i = 1:numel(fields)
+    param.(fields{i}) = param1.(fields{i});  % struct2 values will overwrite struct1 values if overlap occurs
+end
+
+end
+
diff --git a/src/base/solvers/param/parseSolverParameters.m b/src/base/solvers/param/parseSolverParameters.m
index 8198ec37ca..930ea6b602 100644
--- a/src/base/solvers/param/parseSolverParameters.m
+++ b/src/base/solvers/param/parseSolverParameters.m
@@ -101,37 +101,4 @@
         % otherwise use the default parameter
         param.(defaultParams{i,1}) = defaultParams{i,2};
     end
-end
-
-%move following set of parameters from solverOnlyParams to param
-if isfield(solverOnlyParams,'maxConc')
-    param.maxConc = solverOnlyParams.maxConc;
-    solverOnlyParams = rmfield(solverOnlyParams,'maxConc');
-end
-if isfield(solverOnlyParams,'method')
-    param.method = solverOnlyParams.method;
-    solverOnlyParams = rmfield(solverOnlyParams,'method');
-end
-if isfield(solverOnlyParams,'maxUnidirectionalFlux')
-    param.maxUnidirectionalFlux = solverOnlyParams.maxUnidirectionalFlux;
-    solverOnlyParams = rmfield(solverOnlyParams,'maxUnidirectionalFlux');
-end
-if isfield(solverOnlyParams,'minUnidirectionalFlux')
-    param.minUnidirectionalFlux = solverOnlyParams.minUnidirectionalFlux;
-    solverOnlyParams = rmfield(solverOnlyParams,'minUnidirectionalFlux');
-end
-if isfield(solverOnlyParams,'internalNetFluxBounds')
-    param.internalNetFluxBounds = solverOnlyParams.internalNetFluxBounds;
-    solverOnlyParams = rmfield(solverOnlyParams,'internalNetFluxBounds');
-end
-if isfield(solverOnlyParams,'externalNetFluxBounds')
-    param.externalNetFluxBounds = solverOnlyParams.externalNetFluxBounds;
-    solverOnlyParams = rmfield(solverOnlyParams,'externalNetFluxBounds');
-end
-if isfield(solverOnlyParams,'rounding')
-    param.rounding = solverOnlyParams.rounding;
-    solverOnlyParams = rmfield(solverOnlyParams,'rounding');
-end
-if isfield(param,'printLevel')
-    solverOnlyParams.printLevel = param.printLevel -1;
-end
+end
\ No newline at end of file
diff --git a/src/base/solvers/solveCobraLP.m b/src/base/solvers/solveCobraLP.m
index 7da85a4155..35667aa406 100644
--- a/src/base/solvers/solveCobraLP.m
+++ b/src/base/solvers/solveCobraLP.m
@@ -61,7 +61,7 @@
 %                     * .rcost:        Reduced costs, dual solution to :math:`lb <= v <= ub`
 %                     * .dual:         dual solution to `A*v ('E' | 'G' | 'L') b`
 %                     * .solver:       Solver used to solve LP problem
-%                     * .method:    Algorithm used by solver to solve LP problem
+%                     * .lpmethod:    Algorithm used by solver to solve LP problem
 %                     * .stat:         Solver status in standardized form
 %
 %                       * 0 - Infeasible problem
@@ -144,7 +144,7 @@
 end
 
 % support for lifting of ill-scaled models
-if problemTypeParams.lifting == 1
+if isfield(solverParams,'lifting') && solverParams.lifting == 1
     largeNb = 1e4;  % suitable for double precision solvers
     [LPproblem] = reformulate(LPproblem, largeNb, printLevel);
 end
@@ -180,6 +180,9 @@
 if ~isfield(LPproblem, 'modelID')
     LPproblem.modelID = 'aModelID';
 end
+if ~isfield(LPproblem, 'names')
+    LPproblem.names = [];
+end
 
 %too time consuming
 % if any(~isfinite(LPproblem.A),'all')
@@ -197,7 +200,7 @@
 
             
 % extract the problem from the structure
-[A, b, c, lb, ub, csense, osense, modelID] = deal(sparse(LPproblem.A), LPproblem.b, LPproblem.c, LPproblem.lb, LPproblem.ub, LPproblem.csense, LPproblem.osense, LPproblem.modelID);
+[A, b, c, lb, ub, csense, osense, modelID, names] = deal(sparse(LPproblem.A), LPproblem.b, LPproblem.c, LPproblem.lb, LPproblem.ub, LPproblem.csense, LPproblem.osense, LPproblem.modelID,LPproblem.names);
 
 if isfield(LPproblem,'basis') && ~isempty(LPproblem.basis)
     basis = LPproblem.basis;
@@ -222,7 +225,7 @@
 stat = 0;
 origStat = [];
 origStatText = [];
-method = '';
+lpmethod = '';
 
 t_start = clock;
 if isempty(solver)
@@ -545,6 +548,12 @@
 
     case 'glpk'
         %% GLPK
+        %           msglev (default: 1)
+        %                  Level of messages output by solver routines:
+        %                   0 - No output.
+        %                   1 - Error messages only.
+        %                   2 - Normal output.
+        %                   3 - Full output (includes informational messages).
         param.msglev = problemTypeParams.printLevel;  % level of verbosity
         param.tolbnd = problemTypeParams.feasTol;  % tolerance
         param.toldj = problemTypeParams.optTol;  % tolerance
@@ -630,55 +639,9 @@
         w = [];
     case 'mosek'
         % mosek
-        % use msklpopt with full control over all mosek parameters
-        % http://docs.mosek.com/7.0/toolbox/Parameters.html
-        % see also
-        % http://docs.mosek.com/7.0/toolbox/A_guided_tour.html#SEC:VIEWSETPARAM
-        % e.g.
-        % http://docs.mosek.com/7.0/toolbox/MSK_IPAR_OPTIMIZER.html
-
-        %[rcode,res]         = mosekopt('param echo(0)',[],solverParams);
-        
-        % Remove outer function specific parameters to avoid crashing solver interfaces
-        solverParams = mosekParamStrip(solverParams);
-
-        param = solverParams;
-        % only set the print level if not already set via solverParams structure
-        if ~isfield(solverParams, 'MSK_IPAR_LOG')
-            switch problemTypeParams.printLevel
-                case 0
-                    echolev = 0;
-                case 1
-                    echolev = 3;
-                case 2
-                    solverParams.MSK_IPAR_LOG_INTPNT = 1;
-                    solverParams.MSK_IPAR_LOG_SIM = 1;
-                    echolev = 3;
-                otherwise
-                    echolev = 0;
-            end
-            if echolev == 0
-                solverParams.MSK_IPAR_LOG = 0;
-                cmd = ['minimize echo(' int2str(echolev) ')'];
-            else
-                cmd = 'minimize';
-            end
-        end
-
-        %https://docs.mosek.com/8.1/toolbox/solving-linear.html
-        if ~isfield(solverParams, 'MSK_DPAR_INTPNT_TOL_PFEAS')
-            solverParams.MSK_DPAR_INTPNT_TOL_PFEAS=problemTypeParams.feasTol;
-        end
-        if ~isfield(solverParams, 'MSK_DPAR_INTPNT_TOL_DFEAS.')
-            solverParams.MSK_DPAR_INTPNT_TOL_DFEAS=problemTypeParams.feasTol;
-        end
-        %If the feasibility tolerance is changed by the solverParams
-        %struct, this needs to be forwarded to the cobra Params for the
-        %final consistency test!
-        if isfield(solverParams,'MSK_DPAR_INTPNT_TOL_PFEAS')
-            problemTypeParams.feasTol = solverParams.MSK_DPAR_INTPNT_TOL_PFEAS;
-        end
-
+        param = mergeCobraParams(solverParams,problemTypeParams);
+        [cmd,mosekParam] = setMosekParam(param);
+     
         % basis reuse - TODO
         % http://docs.mosek.com/7.0/toolbox/A_guided_tour.html#section-node-_A%20guided%20tour_Advanced%20start%20%28hot-start%29
 
@@ -708,13 +671,12 @@
         %              blx and bux. Note -inf is allowed in blc and blx.
         %              Similarly, inf is allowed in buc and bux.
 
-        
+
         blc = b;
         buc = b;
         buc(csense == 'G') = inf;
         blc(csense == 'L') = -inf;
         
-
         prob.c = osense * c;
         prob.a = A;
         prob.blc     = blc;
@@ -729,8 +691,28 @@
             prob.sol.bas.xx   = basis.xx;
         end
 
-        [rcode,res] = mosekopt(cmd,prob,solverParams);
+        if isfield(param,'saveProb') && param.saveProb
+            formattedTime = datestr(now, 'yyyymmddHHMMSS');
+            LP.cmd=cmd;
+            LP.prob=prob;
+            LP.param=mosekParam;
+            save([formattedTime '_LP_probBeforeMosekopt'],"LP");
+        end
 
+        if isfield(param,'debug') && param.debug==1
+            prob.names = names;
+        end
+
+        if isfield(param,'defaultSolverParam') && param.defaultSolverParam
+            [rcode,res] = mosekopt(cmd,prob);
+        else
+            [rcode,res] = mosekopt(cmd,prob,mosekParam);
+        end
+        if isfield(param,'lpmethod')
+            lpmethod = param.lpmethod;
+        else
+            lpmethod = 'FREE';
+        end
 
         if rcode~=0
             % MSK_RES_TRM_STALL
@@ -745,7 +727,7 @@
         end
             
         %parse mosek result structure
-        [stat,origStat,x,y,yl,yu,z,zl,zu,k,basis,pobjval,dobjval] = parseMskResult(res,solverParams,problemTypeParams.printLevel);
+        [stat,origStat,x,y,yl,yu,z,zl,zu,k,basis,pobjval,dobjval] = parseMskResult(res);
         
         if stat ==1 || stat ==3
             f = c'*x;
@@ -762,10 +744,7 @@
             f = NaN;
             s = NaN*ones(size(A,1),1);
         end
-            
-        if isfield(param,'MSK_IPAR_OPTIMIZER')
-            method=param.MSK_IPAR_OPTIMIZER;
-        end
+        
     case 'mosek_linprog'
         %% mosek
         % if mosek is installed, and the paths are added ahead of matlab's
@@ -822,88 +801,8 @@
     case 'gurobi'
         % Free academic licenses for the Gurobi solver can be obtained from
         % http://www.gurobi.com/html/academic.html
-        % resultgurobi = struct('x',[],'objval',[],'pi',[]);
-
-        %  The params struct contains Gurobi parameters. A full list may be
-        %  found on the Parameter page of the reference manual:
-        %  https://www.gurobi.com/documentation/current/refman/parameter_descriptions.html
-        % MATLAB Parameter Examples
-        % In the MATLAB interface, parameters are passed to Gurobi through a struct. 
-        % To modify a parameter, you create a field in the struct with the appropriate name, 
-        % and set it to the desired value. For example, to set the TimeLimit parameter to 100 you'd do:
-        % 
-        % params.timelimit = 100;
-        % The case of the parameter name is ignored, as are underscores. Thus, you could also do:
-        % params.timeLimit = 100;
-        % ...or...
-        % params.TIME_LIMIT = 100;
-        % All desired parameter changes should be stored in a single struct, which is passed as the second parameter to the gurobi function.
-        param=solverParams;
-
-        % https://www.gurobi.com/documentation/current/refman/method.html
-        % params.method gives the method used to solve continuous models
-        % -1=automatic,
-        %  0=primal simplex,
-        %  1=dual simplex,
-        %  2=barrier,
-        %  3=concurrent,
-        %  4=deterministic concurrent
-        % i.e. params.method     = 1;          % use dual simplex method
-        if isfield(param,'lpmethod')
-            %gurobiAlgorithms = {'AUTOMATIC','PRIMAL','DUAL','BARRIER','CONCURRENT','CONCURRENT_DETERMINISTIC'};
-            % -1=automatic,
-            % 0=primal simplex,
-            % 1=dual simplex,
-            % 2=barrier,
-            % 3=concurrent,
-            % 4=deterministic concurrent
-            switch param.lpmethod
-                case 'AUTOMATIC' 
-                    param.method = -1;
-                case 'PRIMAL'
-                    param.method = 0;
-                case 'DUAL'
-                    param.method = 2;
-                case 'BARRIER'
-                    param.method = 2;
-                case 'CONCURRENT'
-                    param.method = 3;
-                case 'DETERMINISTIC_CONCURRENT'
-                    param.method = 4;
-                otherwise
-                    error('Unrecognised param.lpmethod for gurobi')
-            end
-            param = rmfield(param,'lpmethod');
-        end
-
-        if ~isfield(param,'OutputFlag')
-            switch problemTypeParams.printLevel
-                case 0
-                    param.OutputFlag = 0;
-                    param.DisplayInterval = 1;
-                case 1
-                    param.OutputFlag = 0;
-                    param.DisplayInterval = 1;
-                otherwise
-                    % silent
-                    param.OutputFlag = 0;
-                    param.DisplayInterval = 1;
-            end
-        end
-
-        if isfield(param,'FeasibilityTol')
-            % update tolerance according to actual setting
-            problemTypeParams.feasTol = param.FeasibilityTol;
-        else
-            param.FeasibilityTol = problemTypeParams.feasTol;
-        end
-            
-        if isfield(param,'OptimalityTol')
-            % update tolerance according to actual setting
-            problemTypeParams.optTol = param.OptimalityTol;
-        else      
-            param.OptimalityTol = problemTypeParams.optTol;
-        end
+        param = mergeCobraParams(solverParams,problemTypeParams);
+        gurobiParam = setGurobiParam(param);
 
         gurobiLP.sense(1:length(b),1) = '=';
         gurobiLP.sense(csense == 'L') = '<';
@@ -930,24 +829,14 @@
         %gurobi wants a dense double vector as an objective
         gurobiLP.obj = double(c)+0;%full
 
-
         % basis reuse - Ronan
         if ~isempty(basis)
             gurobiLP.cbasis = full(basis.cbasis);
             gurobiLP.vbasis = full(basis.vbasis);
         end
-        % set the solver specific parameters
-        param = updateStructData(param,solverParams);
 
-%         LPproblem = rmfield(LPproblem,'c');
-%         LPproblem = rmfield(LPproblem,'b');
-%         LPproblem = rmfield(LPproblem,'lb');
-%         LPproblem = rmfield(LPproblem,'ub');
-%         LPproblem = rmfield(LPproblem,'osense');
-%         LPproblem = rmfield(LPproblem,'csense');
-        
         % call the solver
-        resultgurobi = gurobi(gurobiLP,param);
+        resultgurobi = gurobi(gurobiLP,gurobiParam);
 
         % see the solvers original status -Ronan
         origStat = resultgurobi.status;
@@ -1204,7 +1093,7 @@
         [solution,LPprob] = solveCobraLPCPLEX(LPproblem,problemTypeParams.printLevel,1,[],[],minNorm);
         solution.basis = LPprob.LPBasis;
         solution.solver = solver;
-        solution.method = method; % dummy
+        solution.lpmethod = lpmethod; % dummy
      %   solution.slack = [];
         if exist([pwd filesep 'clone1.log'],'file')
             delete('clone1.log')
@@ -1317,7 +1206,7 @@
             %this is the dual to the simple ineequality constraints : reduced costs
             w =  lambda.lower - lambda.upper;
             
-            method = output.method;
+            lpmethod = output.method;
             if 0 %debug
                 disp(method)
                 norm(osense * c  - A' * y - w,inf)
@@ -1364,11 +1253,14 @@
             x = CplexLPproblem.Solution.x;
             if 1
                 f = c'*x;
+                w = CplexLPproblem.Solution.reducedcost;
+                y = CplexLPproblem.Solution.dual;
             else
                 f = CplexLPproblem.Solution.objval; %do not use as it gives the opposite sign for maximise
+                w = osense*CplexLPproblem.Solution.reducedcost;
+                y = osense*CplexLPproblem.Solution.dual;
             end
-            w = osense*CplexLPproblem.Solution.reducedcost;
-            y = osense*CplexLPproblem.Solution.dual;
+
             %res1 = A*solution.full + solution.slack - b;
             s = b - A * x; % output the slack variables
             
@@ -1417,6 +1309,8 @@
             if isfield(CplexLPproblem.Solution ,'dual')
                 y = osense*CplexLPproblem.Solution.dual;
             end
+        elseif origStat==101
+            warning('101')
         else
             stat = -1;
         end
@@ -1645,11 +1539,11 @@
 
 if ~strcmp(solver,'cplex_direct') && ~strcmp(solver,'mps')
     % assign solution
-    t = etime(clock, t_start);
+    time = etime(clock, t_start);
     if ~exist('basis','var'), basis=[]; end
     [solution.full, solution.obj, solution.rcost, solution.dual, solution.slack, ...
-     solution.solver, solution.method, solution.stat, solution.origStat, ...
-     solution.origStatText,solution.time,solution.basis] = deal(x,f,w,y,s,solver,method,stat,origStat,origStatText,t,basis);
+     solution.solver, solution.lpmethod, solution.stat, solution.origStat, ...
+     solution.origStatText,solution.time,solution.basis] = deal(x,f,w,y,s,solver,lpmethod,stat,origStat,origStatText,time,basis);
 elseif strcmp(solver,'mps')
     solution = [];
 end
diff --git a/src/base/solvers/solveCobraQP.m b/src/base/solvers/solveCobraQP.m
index 07c959aba5..12de31c3a5 100644
--- a/src/base/solvers/solveCobraQP.m
+++ b/src/base/solvers/solveCobraQP.m
@@ -62,6 +62,7 @@
 %                       * .slack:       slack variable such that :math:`A*x + s = b`
 %                       * .obj:         Objective value
 %                       * .solver:      Solver used to solve QP problem
+%                       * .qpmethod:    Solver qpmethod used to solve QP problem
 %                       * .origStat:    Original status returned by the specific solver
 %                       * .time:        Solve time in seconds
 %                       * .stat:        Solver status in standardized form (see below)
@@ -323,8 +324,9 @@
             stat = 3; % Solution exists, but either scaling problems or not proven to be optimal
         else %(origStat >= 10)
             stat = -1; % No optimal solution found (time or other limits reached, other infeasibility problems)
-        end
-        
+        end  
+        origStat = cplexProblem.Solution.statusstring;
+
         %Update Tolerance According to actual setting
         problemTypeParams.feasTol = cplexProblem.Param.simplex.tolerances.feasibility.Cur;
         problemTypeParams.optTol = cplexProblem.Param.simplex.tolerances.optimality.Cur;
@@ -372,12 +374,10 @@
         %%
     case 'mosek'
         
-        if problemTypeParams.printLevel>0
-            cmd='minimize';
-        else
-            cmd='minimize echo(0)';
-        end
+        param = mergeCobraParams(solverParams,problemTypeParams);
+        [cmd,mosekParam] = setMosekParam(param);
         
+
         %matching bounds and zero diagonal of F at the same time
         bool = lb == ub & diag(F)==0;
         if any(bool)
@@ -389,38 +389,6 @@
             %}
             warning(['There are ' num2str(nnz(bool)) ' variables that have equal lower and upper bounds, and zero on the diagonal of F.'])
         end
-        
-        param = struct();
-        % Set the printLevel, can be overwritten.
-        if ~isfield(param, 'MSK_IPAR_LOG')
-            switch problemTypeParams.printLevel
-                case 0
-                    echolev = 0;
-                case 1
-                    echolev = 3;
-                case 2
-                    param.MSK_IPAR_LOG_INTPNT = 1;
-                    param.MSK_IPAR_LOG_SIM = 1;
-                    echolev = 3;
-                otherwise
-                    echolev = 0;
-            end
-            if echolev == 0
-                param.MSK_IPAR_LOG = 0;
-                cmd = ['minimize echo(' int2str(echolev) ')'];
-            else
-                cmd = 'minimize';
-            end
-        end
-        %remove parameter fields that mosek does not recognise
-        param.MSK_DPAR_INTPNT_QO_TOL_DFEAS = problemTypeParams.optTol;
-        param.MSK_DPAR_INTPNT_QO_TOL_PFEAS = problemTypeParams.feasTol;
-        
-        %Update with solver Specific Parameter struct
-        param = updateStructData(param,solverParams);
-        %problemTypeParams.feasTol = param.MSK_DPAR_INTPNT_NL_TOL_PFEAS;
-        
-        param = mosekParamStrip(param);
 
         blc = b;
         buc = b;
@@ -438,8 +406,22 @@
         %https://docs.mosek.com/latest/toolbox/data-types.html#prob
         [prob.qosubi,prob.qosubj,prob.qoval]=find(F);
 
-        [rcode,res] = mosekopt(cmd,prob,param);
-      
+        if isfield(param,'saveProb') && param.saveProb
+            formattedTime = datestr(now, 'yyyymmddHHMMSS');
+            QP.cmd=cmd;
+            QP.prob=prob;
+            QP.param=mosekParam;
+            save([formattedTime '_QP_probBeforeMosekopt'],"QP");
+        end
+
+        [rcode,res] = mosekopt(cmd,prob,mosekParam);
+
+        if isfield(param,'qpmethod')
+            qpmethod = param.qpmethod;
+        else
+            qpmethod = 'FREE';
+        end
+
         % stat   Solver status
         %           1   Optimal solution found
         %           2   Unbounded solution
@@ -460,7 +442,7 @@
         end
 
         %parse mosek result structure
-        [stat,origStat,x,y,yl,yu,z,zl,zu,k,basis,pobjval,dobjval] = parseMskResult(res,solverParams,problemTypeParams.printLevel);
+        [stat,origStat,x,y,yl,yu,z,zl,zu,k,basis,pobjval,dobjval] = parseMskResult(res);
         
         %debugging
         if problemTypeParams.printLevel>2
@@ -493,8 +475,7 @@
             f = NaN;
             s = NaN*ones(size(A,1),1);
         end
-        
-        
+
     case 'pdco'
         %-----------------------------------------------------------------------
         % pdco.m: Primal-Dual Barrier Method for Convex Objectives (16 Dec 2008)
@@ -633,74 +614,8 @@
         % Free academic licenses for the Gurobi solver can be obtained from
         % http://www.gurobi.com/html/academic.html
 
-        %  The param struct contains Gurobi parameters. A full list may be
-        %  found on the Parameter page of the reference manual:
-        %  https://www.gurobi.com/documentation/current/refman/parameter_descriptions.html
-        % MATLAB Parameter Examples
-        % In the MATLAB interface, parameters are passed to Gurobi through a struct. 
-        % To modify a parameter, you create a field in the struct with the appropriate name, 
-        % and set it to the desired value. For example, to set the TimeLimit parameter to 100 you'd do:
-        % 
-        % param.timelimit = 100;
-        % The case of the parameter name is ignored, as are underscores. Thus, you could also do:
-        % param.timeLimit = 100;
-        % ...or...
-        % param.TIME_LIMIT = 100;
-        % All desired parameter changes should be stored in a single struct, which is passed as the second parameter to the gurobi function.
-        param=solverParams;
-
-        % param.method gives the method used to solve continuous models
-        % -1=automatic,
-        %  0=primal simplex,
-        %  1=dual simplex,
-        %  2=barrier,
-        %  3=concurrent,
-        %  4=deterministic concurrent
-        % i.e. param.method     = 1;          % use dual simplex method
-        if isfield(param,'lpmethod')
-            %gurobiAlgorithms = {'AUTOMATIC','PRIMAL','DUAL','BARRIER','CONCURRENT','CONCURRENT_DETERMINISTIC'};
-            % -1=automatic,
-            % 0=primal simplex,
-            % 1=dual simplex,
-            % 2=barrier,
-            % 3=concurrent,
-            % 4=deterministic concurrent
-            switch param.lpmethod
-                case 'AUTOMATIC' 
-                    param.method = -1;
-                case 'PRIMAL'
-                    param.method = 0;
-                case 'DUAL'
-                    param.method = 2;
-                case 'BARRIER'
-                    param.method = 2;
-                otherwise
-                    %https://www.gurobi.com/documentation/current/refman/method.html
-                    %Concurrent methods aren't available for QP and QCP. 
-                    error('Unrecognised param.lpmethod for gurobi')
-            end
-            param = rmfield(param,'lpmethod');
-        end
-
-        switch problemTypeParams.printLevel
-            case 0
-                param.OutputFlag = 0;
-                param.DisplayInterval = 1;
-            case problemTypeParams.printLevel>1
-                param.OutputFlag = 1;
-                param.DisplayInterval = 5;
-            otherwise
-                param.OutputFlag = 0;
-                param.DisplayInterval = 1;
-        end
-
-        param.FeasibilityTol = problemTypeParams.feasTol;
-        param.OptimalityTol = problemTypeParams.optTol;
-        %Update param struct with Solver Specific parameters
-        param = updateStructData(param,solverParams);
-        
-        %Update feasTol in case it is changed by the solver Parameters
-        problemTypeParams.feasTol = param.FeasibilityTol;
+        param = mergeCobraParams(solverParams,problemTypeParams);
+        gurobiParam = setGurobiParam(param);
         
         gurobiQP.sense(1:length(b),1) = '=';
         gurobiQP.sense(csense == 'L') = '<';
@@ -730,7 +645,7 @@
         end
         
         try
-            resultgurobi = gurobi(gurobiQP,param);
+            resultgurobi = gurobi(gurobiQP,gurobiParam);
         catch ME
             if contains(ME.message,'Gurobi error 10020: Objective Q not PSD (negative diagonal entry)')
                 warning('%s\n','Gurobi cannot solve a QP problem if it is given a diagonal Q with some of those diagonals equal to zero')
@@ -748,7 +663,7 @@
             
             [x,f,y,w,s] = deal(resultgurobi.x,resultgurobi.objval,osense*resultgurobi.pi,osense*resultgurobi.rc,resultgurobi.slack);
             
-            if problemTypeParams.printLevel>2 %|| 1
+            if param.printLevel>2 %|| 1
                 res1 = A*x + s - b;
                 disp('Check A*x + s - b = 0 (feasiblity):');
                 disp(norm(res1,inf))
@@ -784,7 +699,7 @@
             % if the status becomes 'OPTIMAL', it is unbounded, otherwise it is infeasible.
             gurobiQP.obj(:) = 0;
             gurobiQP.F(:,:) = 0;
-            resultgurobi = gurobi(gurobiQP,param);
+            resultgurobi = gurobi(gurobiQP,gurobiParam);
             if strcmp(resultgurobi.status,'OPTIMAL')
                 stat = 2;
             else
@@ -1047,6 +962,11 @@
 solution.slack = s;
 solution.dual = y;
 solution.rcost = w;
+if exist('qpmethod','var')
+    solution.qpmethod = qpmethod;
+else
+    solution.qpmethod = '';
+end
 if any(contains(solver,'cplex'))
     [ExitText,~] = cplexStatus(solution.origStat);
     solution.origStatText = ExitText;
@@ -1068,7 +988,7 @@
                 fprintf('%s\n',['[' solver '] reports ' solution.origStat ' but Primal optimality condition in solveCobraQP not satisfied, residual = ' num2str(tmp1) ', while feasTol = ' num2str(problemTypeParams.feasTol)])
             else
                 if problemTypeParams.printLevel > 0
-                    fprintf(['\n > [' solver '] Primal optimality condition in solveCobraQP satisfied.']);
+                    fprintf(['\n > [' solver '] Primal optimality condition in solveCobraQP satisfied.\n']);
                 end
             end
         end
diff --git a/test/models b/test/models
index 4ec76c8398..e110c00490 160000
--- a/test/models
+++ b/test/models
@@ -1 +1 @@
-Subproject commit 4ec76c8398e1fb56f1ff67682ea850004ef8f85a
+Subproject commit e110c0049066b4f07d76bea9f31d364e8a1ab977
diff --git a/test/verifiedTests/base/testSolvers/createToyModelForLifting.m b/test/verifiedTests/base/testLifting/createToyModelForLifting.m
similarity index 100%
rename from test/verifiedTests/base/testSolvers/createToyModelForLifting.m
rename to test/verifiedTests/base/testLifting/createToyModelForLifting.m
diff --git a/test/verifiedTests/base/testSolvers/testLiftModel.m b/test/verifiedTests/base/testLifting/testLiftModel.m
similarity index 100%
rename from test/verifiedTests/base/testSolvers/testLiftModel.m
rename to test/verifiedTests/base/testLifting/testLiftModel.m
diff --git a/test/verifiedTests/base/testSolvers/runLPvariousSolvers.m b/test/verifiedTests/base/testSolvers/runLPvariousSolvers.m
index 09f6493a77..30d7385c61 100644
--- a/test/verifiedTests/base/testSolvers/runLPvariousSolvers.m
+++ b/test/verifiedTests/base/testSolvers/runLPvariousSolvers.m
@@ -298,10 +298,10 @@
         else
             tmp_rcost = full(solution{i}.rcost(randrcost));
         end
-        fprintf('%3d%15f%15f%15f%15f%20s\t%30s\n', i, solution{i}.time, solution{i}.obj, tmp_dual, tmp_rcost, solution{i}.solver, solution{i}.method)
+        fprintf('%3d%15f%15f%15f%15f%20s\t%30s\n', i, solution{i}.time, solution{i}.obj, tmp_dual, tmp_rcost, solution{i}.solver, solution{i}.lpmethod)
         all_obj(i) = solution{i}.obj;
     else
-        fprintf('%3d%15f%15f%15f%15f%20s\t%30s\n', i, solution{i}.time, solution{i}.obj, NaN, NaN, solution{i}.solver, solution{i}.method)
+        fprintf('%3d%15f%15f%15f%15f%20s\t%30s\n', i, solution{i}.time, solution{i}.obj, NaN, NaN, solution{i}.solver, solution{i}.lpmethod)
         all_obj(i) = 0.0;
     end
 end
diff --git a/test/verifiedTests/base/testSolvers/testSolveCobraLP.m b/test/verifiedTests/base/testSolvers/testSolveCobraLP.m
index 2acfa54d09..e7d16cb9ee 100644
--- a/test/verifiedTests/base/testSolvers/testSolveCobraLP.m
+++ b/test/verifiedTests/base/testSolvers/testSolveCobraLP.m
@@ -168,8 +168,8 @@
 osenseStr = 'max';
 minNorm = 'zero';
 allowLoops = 1;
-param.zeroNormApprox = 'all';
-solution = optimizeCbModel(model, osenseStr,minNorm, allowLoops, param);
+optimizeCbModel_param.zeroNormApprox = 'all';
+solution = optimizeCbModel(model, osenseStr,minNorm, allowLoops, optimizeCbModel_param);
 assert(solution.f0==1)
 %%
 % change the directory
diff --git a/tutorials b/tutorials
index d87873ff35..ead60d4a2c 160000
--- a/tutorials
+++ b/tutorials
@@ -1 +1 @@
-Subproject commit d87873ff3559a2bb34e02d61b3933ef238eb55b7
+Subproject commit ead60d4a2c63f24d702093dd02b7ac868e14babf