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