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. Only used when minFlux==3 (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 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) Rasmus Agren, 2013-07-16
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. Only used when 0019 % minFlux==3 (opt) 0020 % hsSol hot-start solution for the LP solver. This can 0021 % significantly speed up the process if many similar 0022 % optimization problems are solved iteratively. Only used if 0023 % minFlux is 0 or 1 (opt) 0024 % 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 % hsSolOut solution to be used as hot-start solution (see the input 0035 % parameters). Only used if minFlux is 0 or 1 0036 % 0037 % Usage: [solution hsSolOut]=solveLP(model,minFlux,params,hsSol) 0038 % 0039 % Rasmus Agren, 2013-07-16 0040 % 0041 0042 if nargin<2 0043 minFlux=0; 0044 end 0045 if nargin<3 0046 params.relGap=0.4; 0047 end 0048 if nargin<4 0049 hsSol=[]; 0050 end 0051 0052 %Default return values 0053 hsSolOut=[]; 0054 solution.x=[]; 0055 solution.f=[]; 0056 solution.stat=-1; 0057 0058 %Ignore the hot-start if the previous solution wasn't feasible 0059 if isfield(hsSol,'prosta') 0060 if strfind(hsSol.prosta,'INFEASIBLE') 0061 hsSol=[]; 0062 end 0063 end 0064 0065 % Setup the problem to feed to MOSEK. 0066 prob=[]; 0067 prob.c=model.c*-1; 0068 prob.a=model.S; 0069 prob.blc=model.b(:,1); 0070 %If model.b has two column, then they are for lower/upper bound on the RHS 0071 prob.buc=model.b(:,min(size(model.b,2),2)); 0072 prob.blx=model.lb; 0073 prob.bux=model.ub; 0074 0075 %If hot-start should be used 0076 if ~isempty(hsSol) 0077 prob.sol.bas=hsSol; 0078 params.MSK_IPAR_SIM_HOTSTART=1; 0079 end 0080 0081 %Use MSK_OPTIMIZER_FREE_SIMPLEX. This should not be necessary, but I've 0082 %noticed that the interior point solver is not as good at finding feasible 0083 %solutions. 0084 params.MSK_IPAR_OPTIMIZER='MSK_OPTIMIZER_FREE_SIMPLEX'; 0085 [crap,res] = mosekopt('minimize echo(0)',prob,getMILPParams(params)); 0086 0087 %Check if the problem was feasible and that the solution was optimal 0088 [isFeasible isOptimal]=checkSolution(res); 0089 0090 %If the problem was infeasible using hot-start it is often possible to 0091 %re-solve it without hot-start and get a feasible solution 0092 if ~isFeasible && ~isempty(hsSol) 0093 prob.sol=rmfield(prob.sol,'bas'); 0094 [crap,res] = mosekopt('minimize echo(0)',prob,getMILPParams(params)); 0095 [isFeasible isOptimal]=checkSolution(res); 0096 end 0097 0098 %Return without solution if the problem was infeasible 0099 if ~isFeasible 0100 solution.msg='The problem is infeasible'; 0101 return; 0102 end 0103 if ~isOptimal 0104 solution.msg='The problem is feasible, but not necessarily optimal'; 0105 solution.stat=0; 0106 else 0107 %All is well 0108 solution.stat=1; 0109 solution.msg='Optimal solution found'; 0110 end 0111 0112 %Construct the output structure 0113 if isfield(res.sol,'bas') 0114 solution.x=res.sol.bas.xx; 0115 if minFlux<=1; 0116 hsSolOut=res.sol.bas; 0117 end 0118 solution.f=res.sol.bas.pobjval; 0119 else 0120 %Interior-point. This is not used at the moment 0121 solution.x=res.sol.itr.xx; 0122 solution.f=res.sol.itr.pobjval; 0123 end 0124 0125 %If either the square, the number, or the sum of fluxes should be minimized 0126 %then the objective function value should be fixed before another 0127 %optimization. It is not correct to fix the reactions which participate in 0128 %the objective function to their values in solution.x, as there can be 0129 %multiple solutions with the same objective function value. In addition, this 0130 %approach could result in numerical issues when several fluxes are fixed. 0131 %Instead a new "fake metabolite" is added to the problem. This metabolite 0132 %is produced by each reaction with the stoichiometry that reaction has in 0133 %the objective function. The equality constraint of that "fake metabolite" 0134 %is then set to be at least as good as the objective function value. 0135 if minFlux~=0 0136 model.S=[model.S;prob.c']; 0137 model.mets=[model.mets;'TEMP']; 0138 0139 %If the constraint on the objective function value is exact there is a 0140 %larger risk of numerical errors. However for the quadratic fitting 0141 %intervals are not allowed 0142 if minFlux~=2 0143 if size(model.b,2)==1 0144 model.b=[model.b model.b]; 0145 end 0146 if solution.f<0 0147 model.b=[model.b;[-inf solution.f*0.999999]]; 0148 else 0149 model.b=[model.b;[-inf solution.f*1.000001]]; 0150 end 0151 else 0152 model.b=[model.b;ones(1,size(model.b,2))*solution.f]; 0153 end 0154 0155 switch minFlux 0156 %The sum of fluxes should be minimized 0157 case 1 0158 %Convert the model to the irreversible format 0159 revRxns=find(model.rev); 0160 if ~isempty(revRxns) 0161 iModel=convertToIrrev(model); 0162 else 0163 iModel=model; 0164 end 0165 0166 %Minimize all fluxes 0167 iModel.c(:)=-1; 0168 sol=solveLP(iModel); 0169 0170 %Map back to reversible fluxes 0171 if sol.stat>=0 0172 solution.x=sol.x(1:numel(model.c)); 0173 solution.x(revRxns)=solution.x(revRxns)-sol.x(numel(model.c)+1:end); 0174 else 0175 fprintf('WARNING: Could not solve the problem of minimizing the sum of fluxes. Uses output from original problem\n'); 0176 solution.stat=-2; 0177 end 0178 %The square of fluxes should be minimized. This only works when 0179 %there is no interval on the mass balance constraints (model.b is a 0180 %vector) 0181 case 2 0182 % if size(model.b,2)==1 0183 % qsol=solveQP(model,model.rxns,zeros(numel(model.lb),1)); 0184 % %There is a problem that the solver seldom converges totally in this 0185 % %kind of minimization. Print a warning but use the fluxes 0186 % if any(qsol.x) 0187 % solution.x=qsol.x; 0188 % if qsol.stat==-1 0189 % fprintf('WARNING: The quadratic fitting did not converge\n'); 0190 % end 0191 % else 0192 % fprintf('WARNING: Could not solve the problem of minimizing the square of fluxes. Uses output from linear program\n'); 0193 % end 0194 % else 0195 % fprintf('WARNING: Cannot minimize square of fluxes when size(model.b,2)==2. Uses output from linear program\n'); 0196 % end 0197 fprintf('WARNING: Quadratic solver currently not working. Uses output from original problem\n'); 0198 solution.stat=-2; 0199 %The number of fluxes should be minimized 0200 case 3 0201 [qx,I]=getMinNrFluxes(model,model.rxns,params); 0202 qx(~I)=0; 0203 if any(I) 0204 solution.x=qx; 0205 else 0206 fprintf('WARNING: Could not solve the problem of minimizing the number of fluxes. Uses output from linear program\n'); 0207 solution.stat=-2; 0208 end 0209 end 0210 end 0211 end