optimizeCbModel

PURPOSE ^

optimizeCbModel Solve a flux balance analysis problem

SYNOPSIS ^

function FBAsolution = optimizeCbModel(model,osenseStr, minNorm, allowLoops)

DESCRIPTION ^

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)

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SOURCE CODE ^

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

Generated on Thu 21-Jun-2012 15:39:23 by m2html © 2003