optimizeCbModel Solve a flux balance analysis problem Solves LP problems of the form: max/min c'*v subject to S*v = b lb <= v <= ub FBAsolution = optimizeCbModel(model,osenseStr,minNormFlag) INPUT model (the following fields are required - others can be supplied) S Stoichiometric matrix b Right hand side = dx/dt c Objective coefficients lb Lower bounds ub Upper bounds OPTIONAL INPUTS osenseStr Maximize ('max')/minimize ('min') (opt, default = 'max') minNorm {(0), 'one', > 0 , n x 1 vector}, where [m,n]=size(S); 0 Default, normal LP 'one' Minimise the Taxicab Norm using LP. min |v| s.t. S*v = b c'v = f lb <= v <= ub ----- The remaining options work only with a valid QP solver: ----- > 0 Minimises the Euclidean Norm of internal fluxes. Typically 1e-6 works well. min ||v|| s.t. S*v = b c'v = f lb <= v <= ub n x 1 Forms the diagonal of positive definiate matrix F in the quadratic program min 0.5*v'*F*v st. S*v = b c'*v = f lb <= v <= ub allowLoops {0,(1)} If true, then instead of a conventional FBA, the solver will run an MILP version which does not allow loops in the final solution. Default is true. Runs much slower when set to false. See addLoopLawConstraints.m to for more info. OUTPUT FBAsolution f Objective value x Primal y Dual w Reduced costs s Slacks stat Solver status in standardized form 1 Optimal solution 2 Unbounded solution 0 Infeasible -1 No solution reported (timelimit, numerical problem etc)
0001 function FBAsolution = optimizeCbModel(model,osenseStr, minNorm, allowLoops) 0002 %optimizeCbModel Solve a flux balance analysis problem 0003 % 0004 % Solves LP problems of the form: max/min c'*v 0005 % subject to S*v = b 0006 % lb <= v <= ub 0007 % FBAsolution = optimizeCbModel(model,osenseStr,minNormFlag) 0008 % 0009 %INPUT 0010 % model (the following fields are required - others can be supplied) 0011 % S Stoichiometric matrix 0012 % b Right hand side = dx/dt 0013 % c Objective coefficients 0014 % lb Lower bounds 0015 % ub Upper bounds 0016 % 0017 %OPTIONAL INPUTS 0018 % osenseStr Maximize ('max')/minimize ('min') (opt, default = 'max') 0019 % 0020 % minNorm {(0), 'one', > 0 , n x 1 vector}, where [m,n]=size(S); 0021 % 0 Default, normal LP 0022 % 'one' Minimise the Taxicab Norm using LP. 0023 % min |v| 0024 % s.t. S*v = b 0025 % c'v = f 0026 % lb <= v <= ub 0027 % ----- 0028 % The remaining options work only with a valid QP solver: 0029 % ----- 0030 % > 0 Minimises the Euclidean Norm of internal fluxes. 0031 % Typically 1e-6 works well. 0032 % min ||v|| 0033 % s.t. S*v = b 0034 % c'v = f 0035 % lb <= v <= ub 0036 % n x 1 Forms the diagonal of positive definiate 0037 % matrix F in the quadratic program 0038 % min 0.5*v'*F*v 0039 % st. S*v = b 0040 % c'*v = f 0041 % lb <= v <= ub 0042 % 0043 % allowLoops {0,(1)} If true, then instead of a conventional FBA, 0044 % the solver will run an MILP version which does not allow 0045 % loops in the final solution. Default is true. 0046 % Runs much slower when set to false. 0047 % See addLoopLawConstraints.m to for more info. 0048 % 0049 %OUTPUT 0050 % FBAsolution 0051 % f Objective value 0052 % x Primal 0053 % y Dual 0054 % w Reduced costs 0055 % s Slacks 0056 % stat Solver status in standardized form 0057 % 1 Optimal solution 0058 % 2 Unbounded solution 0059 % 0 Infeasible 0060 % -1 No solution reported (timelimit, numerical problem etc) 0061 % 0062 0063 % Markus Herrgard 9/16/03 0064 % Ronan Fleming 4/25/09 Option to minimises the Euclidean Norm of internal 0065 % fluxes using 'cplex_direct' solver 0066 % Ronan Fleming 7/27/09 Return an error if any imputs are NaN 0067 % Ronan Fleming 10/24/09 Fixed 'E' for all equality constraints 0068 % Jan Schellenberger MILP option to remove flux around loops 0069 % Ronan Fleming 12/07/09 Reworked minNorm parameter option to allow 0070 % the full range of approaches for getting 0071 % rid of net flux around loops. 0072 % Jan Schellenberger 2/3/09 fixed bug with .f being set incorrectly 0073 % when minNorm was set. 0074 % Nathan Lewis 12/2/10 Modified code to allow for inequality 0075 % constraints. 0076 % Ronan Fleming 12/03/10 Minor changes to the internal handling of global parameters. 0077 %% Process arguments and set up problem 0078 0079 if exist('osenseStr', 'var') 0080 if isempty(osenseStr) 0081 osenseStr = 'max'; 0082 end 0083 else 0084 osenseStr = 'max'; 0085 end 0086 0087 if exist('minNorm', 'var') 0088 if isempty(minNorm) 0089 minNorm = false; 0090 changeOK = changeCobraSolverParams('LP','minNorm',minNorm); 0091 else 0092 changeOK = changeCobraSolverParams('LP','minNorm',minNorm); 0093 end 0094 else 0095 minNorm = false; 0096 changeOK = changeCobraSolverParams('LP','minNorm',minNorm); 0097 end 0098 if exist('allowLoops', 'var') 0099 if isempty(allowLoops) 0100 allowLoops = true; 0101 end 0102 else 0103 allowLoops = true; 0104 end 0105 0106 [minNorm, printLevel, primalOnlyFlag, saveInput] = getCobraSolverParams('LP',{'minNorm','printLevel','primalOnly','saveInput'}); 0107 0108 0109 % if exist('minNorm', 'var') 0110 % if isempty(minNorm) 0111 % minNorm = false; 0112 % end 0113 % else 0114 % minNorm = false; 0115 % end 0116 % if exist('allowLoops', 'var') 0117 % if isempty(allowLoops) 0118 % allowLoops = true; 0119 % end 0120 % else 0121 % allowLoops = true; 0122 % end 0123 % 0124 % 0125 % global CBT_LP_PARAMS 0126 % if (exist('CBT_LP_PARAMS', 'var')) 0127 % if isfield(CBT_LP_PARAMS, 'objTol') 0128 % tol = CBT_LP_PARAMS.objTol; 0129 % else 0130 % tol = 1e-6; 0131 % end 0132 % if isfield(CBT_LP_PARAMS, 'primalOnly') 0133 % primalOnlyFlag = CBT_LP_PARAMS.primalOnly; 0134 % else 0135 % primalOnlyFlag = false; 0136 % end 0137 % if isfield(CBT_LP_PARAMS, 'printLevel') 0138 % printLevel = CBT_LP_PARAMS.printLevel; 0139 % else 0140 % printLevel = 0; 0141 % end 0142 % else 0143 % tol = 1e-6; 0144 % primalOnlyFlag = false; 0145 % printLevel = 0; 0146 % end 0147 0148 % Figure out objective sense 0149 if strcmpi(osenseStr,'max') 0150 LPproblem.osense = -1; 0151 else 0152 LPproblem.osense = +1; 0153 end 0154 0155 % this is dangerous... if model does not have S, it should not be called in 0156 % this function. 0157 % if ~isfield(model,'S') 0158 % model.S=model.A; 0159 % end 0160 0161 [nMets,nRxns] = size(model.S); 0162 0163 % add csense 0164 %Doing this makes csense a double array. Totally smart design move. 0165 %LPproblem.csense = []; 0166 if ~isfield(model,'csense') 0167 % If csense is not declared in the model, assume that all 0168 % constraints are equalities. 0169 LPproblem.csense(1:nMets,1) = 'E'; 0170 else % if csense is in the model, move it to the lp problem structure 0171 if length(model.csense)~=nMets, 0172 warning('Length of csense is invalid! Defaulting to equality constraints.') 0173 LPproblem.csense(1:nMets,1) = 'E'; 0174 else 0175 model.csense = columnVector(model.csense); 0176 LPproblem.csense = model.csense; 0177 end 0178 end 0179 0180 0181 % Fill in the RHS vector if not provided 0182 if (~isfield(model,'b')) 0183 LPproblem.b = zeros(size(model.S,1),1); 0184 else 0185 LPproblem.b = model.b; 0186 end 0187 0188 % Rest of the LP problem 0189 LPproblem.A = model.S; 0190 LPproblem.c = model.c; 0191 LPproblem.lb = model.lb; 0192 LPproblem.ub = model.ub; 0193 0194 %Double check that all inputs are valid: 0195 if ~(verifyCobraProblem(LPproblem, [], [], false) == 1) 0196 warning('invalid problem'); 0197 return; 0198 end 0199 0200 %% 0201 t1 = clock; 0202 % Solve initial LP 0203 if allowLoops 0204 solution = solveCobraLP(LPproblem); 0205 else 0206 MILPproblem = addLoopLawConstraints(LPproblem, model, 1:nRxns); 0207 solution = solveCobraMILP(MILPproblem); 0208 end 0209 0210 if (solution.stat ~= 1) % check if initial solution was successful. 0211 if printLevel>0 0212 warning('Optimal solution was not found'); 0213 end 0214 FBAsolution.f = 0; 0215 FBAsolution.x = []; 0216 FBAsolution.stat = solution.stat; 0217 FBAsolution.origStat = solution.origStat; 0218 FBAsolution.solver = solution.solver; 0219 FBAsolution.time = etime(clock, t1); 0220 return; 0221 end 0222 0223 objective = solution.obj; % save for later use. 0224 0225 if strcmp(minNorm, 'one') 0226 % Minimize the absolute value of fluxes to 'avoid' loopy solutions 0227 % Solve secondary LP to minimize one-norm of |v| 0228 % Set up the optimization problem 0229 % min sum(delta+ + delta-) 0230 % 1: S*v1 = 0 0231 % 3: delta+ >= -v1 0232 % 4: delta- >= v1 0233 % 5: c'v1 >= f (optimal value of objective) 0234 % 0235 % delta+,delta- >= 0 0236 LPproblem2.A = [model.S sparse(nMets,2*nRxns); 0237 speye(nRxns,nRxns) speye(nRxns,nRxns) sparse(nRxns,nRxns); 0238 -speye(nRxns,nRxns) sparse(nRxns,nRxns) speye(nRxns,nRxns); 0239 model.c' sparse(1,2*nRxns)]; 0240 LPproblem2.c = [zeros(nRxns,1);ones(2*nRxns,1)]; 0241 LPproblem2.lb = [model.lb;zeros(2*nRxns,1)]; 0242 LPproblem2.ub = [model.ub;10000*ones(2*nRxns,1)]; 0243 LPproblem2.b = [LPproblem.b;zeros(2*nRxns,1);solution.obj]; 0244 if ~isfield(model,'csense') 0245 % If csense is not declared in the model, assume that all 0246 % constraints are equalities. 0247 LPproblem2.csense(1:nMets) = 'E'; 0248 else % if csense is in the model, move it to the lp problem structure 0249 if length(model.csense)~=nMets, 0250 warning('Length of csense is invalid! Defaulting to equality constraints.') 0251 LPproblem2.csense(1:nMets) = 'E'; 0252 else 0253 LPproblem2.csense = columnVector(model.csense); 0254 end 0255 end 0256 LPproblem2.csense((nMets+1):(nMets+2*nRxns)) = 'G'; 0257 LPproblem2.csense(nMets+2*nRxns+1) = 'G'; 0258 LPproblem2.csense = columnVector(LPproblem2.csense); 0259 LPproblem2.osense = 1; 0260 % Re-solve the problem 0261 if allowLoops 0262 solution = solveCobraLP(LPproblem2); %,printLevel,minNorm); 0263 solution.dual = []; % slacks and duals will not be valid for this computation. 0264 solution.rcost = []; 0265 else 0266 MILPproblem2 = addLoopLawConstraints(LPproblem, model, 1:nRxns); 0267 solution = solveCobraMILP(MILPproblem2); 0268 end 0269 elseif length(minNorm)> 1 || minNorm > 0 0270 % quadratic minimization of the norm. 0271 % set previous optimum as constraint. 0272 LPproblem.A = [LPproblem.A; 0273 LPproblem.c']; 0274 LPproblem.csense(end+1) = 'E'; 0275 if nnz(LPproblem.c)>1 0276 error('Code assumes only one non-negative coefficient in linear part of objective'); 0277 end 0278 LPproblem.b = [LPproblem.b;solution.full(LPproblem.c~=0)]; 0279 LPproblem.c = zeros(size(LPproblem.c)); % no need for c anymore. 0280 %Minimise Euclidean norm using quadratic programming 0281 if length(minNorm)==1 0282 minNorm=ones(nRxns,1)*minNorm; 0283 end 0284 LPproblem.F = spdiags(minNorm,0,nRxns,nRxns); 0285 %quadratic optimization 0286 if allowLoops 0287 solution = solveCobraQP(LPproblem); 0288 else 0289 MIQPproblem = addLoopLawConstraints(LPproblem, model, 1:nRxns); 0290 solution = solveCobraMIQP(MIQPproblem); 0291 end 0292 % if isempty(solution.full) 0293 % % QP problem did not work. This will return empty structure later. 0294 % else 0295 % %dont include dual variable to additional constraint 0296 % %solution.dual=solution.dual(1:end-1,1); 0297 % end 0298 end 0299 0300 % Store results 0301 if (solution.stat == 1) 0302 %solution found. 0303 FBAsolution.x = solution.full(1:nRxns); 0304 0305 %this line IS necessary. 0306 FBAsolution.f = model.c'*solution.full(1:nRxns); %objective from original optimization problem. 0307 if abs(FBAsolution.f - objective) > .01 0308 display('warning: objective appears to have changed while performing secondary optimization (minNorm)'); 0309 end 0310 0311 if (~primalOnlyFlag && allowLoops && any(~minNorm)) % rcost/dual only correct if not doing minNorm 0312 FBAsolution.y = solution.dual; 0313 FBAsolution.w = solution.rcost; 0314 end 0315 else 0316 %some sort of error occured. 0317 if printLevel>0 0318 warning('Optimal solution was not found'); 0319 end 0320 FBAsolution.f = 0; 0321 FBAsolution.x = []; 0322 end 0323 0324 FBAsolution.stat = solution.stat; 0325 FBAsolution.origStat = solution.origStat; 0326 FBAsolution.solver = solution.solver; 0327 FBAsolution.time = etime(clock, t1); 0328