solveLP
Solves a linear programming problem
model a model structure
minFlux determines if a second optimization should be performed
in order to get rid of loops in the flux distribution
0: no such optimization is performed
1: the sum of abs(fluxes) is minimized. This is the
fastest way of getting rid of loops
2: the square of fluxes is minimized. This tends to
distribute fluxes across iso-enzymes, which results in a
larger number of reactions being used
3: the number of fluxes is minimized. This can result
in the flux distributions that are the easiest to
interpret. Note that this optimization can be very slow
(opt, default 0)
params parameter structure as used by getMILPParams (opt)
hsSol hot-start solution for the LP solver. This can
significantly speed up the process if many similar
optimization problems are solved iteratively. Only used if
minFlux is 0 or 1 (opt)
solution
f objective value
x primal (flux distribution)
stat exit flag
1: the optimization terminated successfully
0: the solution is feasible, but not necessarily optimal
-1: no feasible solution found
-2: solution found, but flux minimization failed
msg string describing the status of the optimization
sPrice shadow price (only reported if minFlux was 0)
rCost reduced cost (only reported if minFlux was 0)
hsSolOut solution to be used as hot-start solution (see the input
parameters). Only used if minFlux is 0 or 1
Usage: [solution, hsSolOut]=solveLP(model,minFlux,params,hsSol)0001 function [solution, hsSolOut]=solveLP(model,minFlux,params,hsSol) 0002 % solveLP 0003 % Solves a linear programming problem 0004 % 0005 % model a model structure 0006 % minFlux determines if a second optimization should be performed 0007 % in order to get rid of loops in the flux distribution 0008 % 0: no such optimization is performed 0009 % 1: the sum of abs(fluxes) is minimized. This is the 0010 % fastest way of getting rid of loops 0011 % 2: the square of fluxes is minimized. This tends to 0012 % distribute fluxes across iso-enzymes, which results in a 0013 % larger number of reactions being used 0014 % 3: the number of fluxes is minimized. This can result 0015 % in the flux distributions that are the easiest to 0016 % interpret. Note that this optimization can be very slow 0017 % (opt, default 0) 0018 % params parameter structure as used by getMILPParams (opt) 0019 % hsSol hot-start solution for the LP solver. This can 0020 % significantly speed up the process if many similar 0021 % optimization problems are solved iteratively. Only used if 0022 % minFlux is 0 or 1 (opt) 0023 % 0024 % solution 0025 % f objective value 0026 % x primal (flux distribution) 0027 % stat exit flag 0028 % 1: the optimization terminated successfully 0029 % 0: the solution is feasible, but not necessarily optimal 0030 % -1: no feasible solution found 0031 % -2: solution found, but flux minimization failed 0032 % msg string describing the status of the optimization 0033 % sPrice shadow price (only reported if minFlux was 0) 0034 % rCost reduced cost (only reported if minFlux was 0) 0035 % hsSolOut solution to be used as hot-start solution (see the input 0036 % parameters). Only used if minFlux is 0 or 1 0037 % 0038 % Usage: [solution, hsSolOut]=solveLP(model,minFlux,params,hsSol) 0039 0040 if nargin<2 0041 minFlux=0; 0042 end 0043 if nargin<3 0044 params.relGap=0.4; 0045 end 0046 if nargin<4 0047 hsSol=[]; 0048 end 0049 0050 %Default return values 0051 hsSolOut=[]; 0052 solution.x=[]; 0053 solution.f=[]; 0054 solution.stat=-1; 0055 0056 %Check for valid lb and ub 0057 if ~isnumeric([model.lb,model.ub]) 0058 invalidBound = true(numel(model.lb),1); 0059 else 0060 invalidBound = false(numel(model.lb),1); 0061 end 0062 invalidBound = invalidBound | model.lb>model.ub; 0063 invalidBound = invalidBound | logical(sum(isnan([model.lb,model.ub]),2)); 0064 if any(invalidBound) 0065 error(['Invalid set of bounds for reaction(s): ', strjoin(model.rxns(invalidBound),', '), '.']) 0066 end 0067 %Check for valid objective 0068 if ~isnumeric(model.c) || any(isnan(model.c)) 0069 error('Invalid defintion of objective function in model.c.') 0070 end 0071 %Check for valid S-matrix 0072 invalidS = ~isfinite(model.S); 0073 if any(any(invalidS)) 0074 error(['Invalid coefficients defined for reaction(s): ', strjoin(model.rxns(any(invalidS)),', '), '.']) 0075 end 0076 0077 %Ignore the hot-start if the previous solution wasn't feasible 0078 if isfield(hsSol,'stat') 0079 if hsSol.stat<1 0080 hsSol=[]; 0081 end 0082 end 0083 0084 % Setup the linear problem for glpk 0085 prob = []; 0086 prob.c = [model.c*-1;zeros(size(model.S,1),1)]; 0087 prob.b = zeros(size(model.S,1), 1); 0088 prob.lb = [model.lb; model.b(:,1)]; 0089 prob.ub = [model.ub; model.b(:,min(size(model.b,2),2))]; 0090 prob.csense = repmat('E', 1, length(prob.b)); % Equality constraints 0091 prob.osense = 1; 0092 prob.vartype = repmat('C', 1, length(prob.c)); % Fluxes are continuous values 0093 prob.A = [model.S -speye(size(model.S,1))]; 0094 prob.a = model.S; 0095 0096 %If hot-start should be used 0097 if ~isempty(hsSol) 0098 prob.vbasis=hsSol.vbasis; 0099 prob.cbasis=hsSol.cbasis; 0100 end 0101 0102 % Parse the problem to the LP solver 0103 res = optimizeProb(prob,params); 0104 0105 %Check if the problem was feasible and that the solution was optimal 0106 [isFeasible, isOptimal]=checkSolution(res); 0107 0108 %If the problem was infeasible using hot-start it is often possible to 0109 %re-solve it without hot-start and get a feasible solution 0110 if ~isFeasible && ~isempty(hsSol) 0111 prob=rmfield(prob,{'vbasis','cbasis'}); 0112 res=optimizeProb(prob,params); 0113 [isFeasible, isOptimal]=checkSolution(res); 0114 end 0115 0116 %Return without solution if the problem was infeasible 0117 if ~isFeasible 0118 solution.msg='The problem is infeasible'; 0119 return; 0120 end 0121 if ~isOptimal 0122 solution.msg='The problem is feasible, but not necessarily optimal'; 0123 solution.stat=0; 0124 else 0125 %All is well 0126 solution.stat=1; 0127 solution.msg='Optimal solution found'; 0128 end 0129 0130 %Construct the output structure 0131 if isfield(res,'full') 0132 solution.x=res.full; 0133 if minFlux<=1 0134 if(isfield(res,'vbasis')) % gurobi uses vbasis and cbasis as hotstart 0135 hsSolOut.vbasis=res.vbasis; 0136 hsSolOut.cbasis=res.cbasis; 0137 end 0138 end 0139 solution.f=res.obj; 0140 if isfield(res,'dual') 0141 solution.sPrice=res.dual; 0142 else 0143 solution.sPrice=[]; 0144 end 0145 if isfield(res,'rcost') 0146 solution.rCost=res.rcost(1:numel(model.rxns)); 0147 else 0148 solution.rCost=[]; 0149 end 0150 0151 %If either the square, the number, or the sum of fluxes should be minimized 0152 %then the objective function value should be fixed before another 0153 %optimization. It is not correct to fix the reactions which participate in 0154 %the objective function to their values in solution.x, as there can be 0155 %multiple solutions with the same objective function value. In addition, 0156 %this approach could result in numerical issues when several fluxes are 0157 %fixed. Instead a new "fake metabolite" is added to the problem. This 0158 %metabolite is produced by each reaction with the stoichiometry that 0159 %reaction has in the objective function. The equality constraint of that 0160 %"fake metabolite" is then set to be at least as good as the objective 0161 %function value. 0162 if minFlux~=0 0163 model.S=[model.S;(model.c*-1)']; 0164 model.mets=[model.mets;'TEMP']; 0165 0166 %If the constraint on the objective function value is exact there is a 0167 %larger risk of numerical errors. However for the quadratic fitting 0168 %intervals are not allowed 0169 if minFlux~=2 0170 if size(model.b,2)==1 0171 model.b=[model.b model.b]; 0172 end 0173 if solution.f<0 0174 model.b=[model.b;[-inf solution.f*0.999999]]; 0175 else 0176 model.b=[model.b;[-inf solution.f*1.000001]]; 0177 end 0178 else 0179 model.b=[model.b;ones(1,size(model.b,2))*solution.f]; 0180 end 0181 0182 switch minFlux 0183 %The sum of fluxes should be minimized 0184 case 1 0185 %Convert the model to the irreversible format 0186 revRxns=find(model.rev); 0187 if ~isempty(revRxns) 0188 iModel=convertToIrrev(model); 0189 else 0190 iModel=model; 0191 end 0192 %Minimize all fluxes 0193 iModel.c(:)=-1; 0194 sol=solveLP(iModel,0,params); 0195 0196 %Map back to reversible fluxes 0197 if sol.stat>=0 0198 solution.x=sol.x(1:numel(model.c)); 0199 solution.x(revRxns)=solution.x(revRxns)-sol.x(numel(model.c)+1:end); 0200 else 0201 EM='Could not solve the problem of minimizing the sum of fluxes. Uses output from original problem'; 0202 dispEM(EM,false); 0203 solution.stat=-2; 0204 end 0205 0206 case 2 0207 error('%Minimization of square of fluxes is deprecated') 0208 0209 %The number of fluxes should be minimized 0210 case 3 0211 [qx,I]=getMinNrFluxes(model,model.rxns,params); 0212 qx(~I)=0; 0213 if any(I) 0214 solution.x=qx; 0215 else 0216 EM='Could not solve the problem of minimizing the number of fluxes. Uses output from linear program'; 0217 dispEM(EM,false); 0218 solution.stat=-2; 0219 end 0220 end 0221 solution=rmfield(solution,{'sPrice','rCost'}); 0222 end 0223 end