Home > RAVEN > solveLP.m

solveLP

PURPOSE ^

solveLP

SYNOPSIS ^

function [solution hsSolOut]=solveLP(model,minFlux,params,hsSol)

DESCRIPTION ^

 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

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SOURCE CODE ^

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

Generated on Tue 16-Jul-2013 21:50:02 by m2html © 2005