createBilevelMILPproblem

PURPOSE ^

createBilevelMILPProblem Creates the necessary matrices and vectors to solve a bilevel MILP with designated inner

SYNOPSIS ^

function bilevelMILPproblem = createBilevelMILPproblem(model,cLinear,cInteger,selRxns,selRxnMatch,constrOpt,measOpt,options,selPrevSol)

DESCRIPTION ^

createBilevelMILPProblem Creates the necessary matrices and vectors to solve a bilevel MILP with designated inner
and outer problem objective functions

 bilevelMILPproblem = createBilevelMILPProblem(model,cLinear,cInteger,selRxns,...
                                               selRxnMatch,constrOpt,measOpt,options,selPrevSol);

INPUTS
 model        Model in irreversible format
 cLinear       Objective for linear part of the MILP problem (i.e. for fluxes)
 cInteger      Objective for integer part of the MILP problem
 selRxns       Reactions that participate in the integer part (e.g. ones
               that can be deleted) (in the form [0 0 1 0 0 1 0 1 1 0 1])
 selRxnMatch  Matching of the forward and reverse parts
 constrOpt    Constraint options (optional)
  rxnInd
  values
  sense
 measOpt      Measured flux options (optional)
  rxnSel
  values
  weights
 options      General options
  vMax         Maximal/minimal value for cont variables
  numDel       Number of deletions
  numDelSense  # of ko <=/=/>= K (L/E/G)
 selPrevSol   Previous solutions (optional)

OUTPUT
 bilevelMILPproblem
   A             LHS matrix 
   b             RHS
   c             Objective
   csense        Constraint types
   lb            Lower bounds
   ub            Upper bounds
   vartype       Variable types
   contSolInd    Allows selecting the continuous solution (i.e. fluxes)
   intsSolInd    Allows selecting the integer solution (i.e. what
                 reactions are used)

 Outputs suitable for feeding to feeding into a MILP solver

 Variables are in order: v(n), y(n_int), u_eq(m), u_u(n), u_z(n_int) u_b(n_ic) s_mm(n_m) s_mp(n_m)

 1/22/07 Add new interface to allow inclusion in the COBRA Toolbox
 
 5/24/05 Fixed problem with reversible measured fluxes - requires changing
 the interpretation of the sel_m input vector (can now have both positive
 & negative entries) MH

 5/27/05 Added ability to rule out particular integer solutions (i.e.
 existing solutions) MH

 6/10/05 Added ability to delimit fluxes to a certain range

 Markus Herrgard 5/27/05

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SOURCE CODE ^

0001 function bilevelMILPproblem = createBilevelMILPproblem(model,cLinear,cInteger,selRxns,...
0002     selRxnMatch,constrOpt,measOpt,options,selPrevSol)
0003 %createBilevelMILPProblem Creates the necessary matrices and vectors to solve a bilevel MILP with designated inner
0004 %and outer problem objective functions
0005 %
0006 % bilevelMILPproblem = createBilevelMILPProblem(model,cLinear,cInteger,selRxns,...
0007 %                                               selRxnMatch,constrOpt,measOpt,options,selPrevSol);
0008 %
0009 %INPUTS
0010 % model        Model in irreversible format
0011 % cLinear       Objective for linear part of the MILP problem (i.e. for fluxes)
0012 % cInteger      Objective for integer part of the MILP problem
0013 % selRxns       Reactions that participate in the integer part (e.g. ones
0014 %               that can be deleted) (in the form [0 0 1 0 0 1 0 1 1 0 1])
0015 % selRxnMatch  Matching of the forward and reverse parts
0016 % constrOpt    Constraint options (optional)
0017 %  rxnInd
0018 %  values
0019 %  sense
0020 % measOpt      Measured flux options (optional)
0021 %  rxnSel
0022 %  values
0023 %  weights
0024 % options      General options
0025 %  vMax         Maximal/minimal value for cont variables
0026 %  numDel       Number of deletions
0027 %  numDelSense  # of ko <=/=/>= K (L/E/G)
0028 % selPrevSol   Previous solutions (optional)
0029 %
0030 %OUTPUT
0031 % bilevelMILPproblem
0032 %   A             LHS matrix
0033 %   b             RHS
0034 %   c             Objective
0035 %   csense        Constraint types
0036 %   lb            Lower bounds
0037 %   ub            Upper bounds
0038 %   vartype       Variable types
0039 %   contSolInd    Allows selecting the continuous solution (i.e. fluxes)
0040 %   intsSolInd    Allows selecting the integer solution (i.e. what
0041 %                 reactions are used)
0042 %
0043 % Outputs suitable for feeding to feeding into a MILP solver
0044 %
0045 % Variables are in order: v(n), y(n_int), u_eq(m), u_u(n), u_z(n_int) u_b(n_ic) s_mm(n_m) s_mp(n_m)
0046 %
0047 % 1/22/07 Add new interface to allow inclusion in the COBRA Toolbox
0048 %
0049 % 5/24/05 Fixed problem with reversible measured fluxes - requires changing
0050 % the interpretation of the sel_m input vector (can now have both positive
0051 % & negative entries) MH
0052 %
0053 % 5/27/05 Added ability to rule out particular integer solutions (i.e.
0054 % existing solutions) MH
0055 %
0056 % 6/10/05 Added ability to delimit fluxes to a certain range
0057 %
0058 % Markus Herrgard 5/27/05
0059 
0060 % Convert inputs
0061 
0062 % clp           Objective for continuous variables in the outer problem
0063 % cip           Objective for integer variables in the outer problem
0064 % S             Stoichiometric matrix (or equality constraint matrix)
0065 % ub            Upper bounds for inner problem
0066 % c             Objective for inner problem
0067 % sel_int       Varibles (inner) corresponding to the integers (in the form [0 0 1 0 0 1 0 0 0 1 1 0 1 1])
0068 % rev_int       Matchings for reversible varibles corresponding to integers (in the form [1 2;4 5])
0069 % ind_ic        Varibles for which there are additional constraints (in the
0070 % form [1 4 4 5 9]);
0071 % b_ic          Additional constraints (in the form [1.5 2.3 3.2 3.3 1.1])
0072 % csense_ic     Directions of additional constraints (in the form 'EGLLG')
0073 % sel_m         Select measured fluxes (in the form [1 0 0 1 -1 0 0 0 0 0 0 0 0])
0074 % b_m           Values for measured fluxes (in the form [2.9 3.0])
0075 % wt_m          Weights for measured fluxes (in the form [0.1 0.9])
0076 % H             Maximal/minimal value for cont variables
0077 % K             Desired number of reaction knockouts
0078 % ksense        # of ko <=/=/>= K (L/E/G)
0079 % sel_int_prev  Varibles corresponding to the previously selected integers
0080 %               (optional). Same form/dimension as sel_int;
0081 
0082 S = model.S;
0083 ub = model.ub;
0084 c = model.c;
0085 clp = cLinear;
0086 cip = cInteger;
0087 sel_int = selRxns;
0088 rev_int = selRxnMatch;
0089 if (~isempty(constrOpt))
0090     ind_ic = constrOpt.rxnInd;
0091     b_ic = constrOpt.values;
0092     csense_ic = constrOpt.sense;
0093 else
0094     ind_ic = [];
0095     b_ic = [];
0096     csense_ic = [];
0097 end
0098 if (~isempty(measOpt))
0099     sel_m = measOpt.rxnSel;
0100     b_m = measOpt.values;
0101     wt_m= measOpt.weights;
0102 else
0103     sel_m = [];
0104     b_m = [];
0105     wt_m = [];
0106 end
0107 H = options.vMax;
0108 K = options.numDel;
0109 ksense = options.numDelSense;
0110 if (nargin < 9)
0111     sel_int_prev = [];
0112 else
0113     sel_int_prev = selPrevSol;
0114 end
0115 
0116 % Dimensions of the primal problem
0117 [m,n] = size(S);
0118 % Number of integer variables
0119 n_int = sum(sel_int);
0120 % Number of inner problem "special" constraints
0121 n_ic = length(ind_ic);
0122 sel_ic = zeros(length(sel_int),1);
0123 sel_ic(ind_ic) = 1;
0124 % Number of inner problem primal varibles with no integers associated with
0125 % them & not part of the special constraint set
0126 n_nint = n - n_int - sum(sel_ic);
0127 % Number of reversible rxns with integer
0128 [n_intr,tmp] = size(rev_int);
0129 % Number of measured fluxes (count only positive directions)
0130 n_m = sum(sel_m == 1);
0131 
0132 % Helper arrays for extracting solutions
0133 sel_cont_sol = [1:n];
0134 sel_int_sol = [n+1:n+n_int];
0135 
0136 % Set variable types
0137 vartype_bl(1:2*n+2*n_int+m+n_ic+2*n_m) = 'C';
0138 vartype_bl(n+1:n+n_int) = 'B';
0139 
0140 % Set upper/lower bounds
0141 lb_bl = [zeros(n+n_int,1); -H*ones(m,1); zeros(n,1); -H*ones(n_int,1)];
0142 ub_bl = [H*ones(2*n+2*n_int+m,1)];
0143 % Handle inner problem constraints
0144 for i = 1:n_ic
0145     % Don't constraint equality constrained ones or
0146     if  (strcmp(csense_ic(i),'E'))
0147         lb_bl = [lb_bl; -H];
0148         ub_bl = [ub_bl; H];
0149     % Lower bound constraints
0150     elseif (strcmp(csense_ic(i),'L'))
0151         lb_bl = [lb_bl; 0];
0152         ub_bl = [ub_bl; H];
0153     % Upper bound constraints
0154     elseif (strcmp(csense_ic(i),'G'))
0155         lb_bl = [lb_bl; -H];
0156         ub_bl = [ub_bl; 0];
0157     else
0158         lb_bl = [lb_bl; -H];
0159         ub_bl = [ub_bl; H];
0160     end
0161 end
0162 % Slacks for measured fluxes
0163 lb_bl = [lb_bl; zeros(2*n_m,1)];
0164 ub_bl = [ub_bl; H*ones(2*n_m,1)];
0165 ub_bl(vartype_bl == 'B') = 1;
0166 
0167 % Set bilevel objective
0168 c_bl = zeros(2*n+2*n_int+m+n_ic+2*n_m,1);
0169 c_bl(1:n) = clp;
0170 c_bl(n+1:n+n_int) = cip;
0171 %c_bl(2*n+2*n_int+m+n_ic+1:2*n+2*n_int+m+n_ic+2*n_m) = 1;
0172 c_bl(2*n+2*n_int+m+n_ic+1:2*n+2*n_int+m+n_ic+n_m) = wt_m;
0173 c_bl(2*n+2*n_int+m+n_ic+n_m+1:2*n+2*n_int+m+n_ic+2*n_m) = wt_m;
0174 
0175 %********************************
0176 % Generate the constraint matrix
0177 %********************************
0178 
0179 % Create necessary integer matrices (for SF/n_int sets)
0180 Iint = selMatrix(sel_int);
0181 Inint = selMatrix(~sel_int & ~sel_ic);
0182 Iic = sparse(n_ic,n);
0183 for i = 1:length(ind_ic)
0184     idx_ic = ind_ic(i);
0185     Iic(i,idx_ic) = 1;
0186 end
0187 Im = selMatrix(sel_m);
0188 
0189 % S*v = 0
0190 A_bl = [S sparse(m,2*n_int+n+m+n_ic+2*n_m)];
0191 b_bl = zeros(m,1);
0192 csense_bl(1:m) = 'E';
0193 
0194 % v_ic >= b_ic
0195 A_bl = [A_bl; Iic sparse(n_ic,2*n_int+m+n+n_ic+2*n_m)];
0196 b_bl = [b_bl; b_ic];
0197 csense_bl(end+1:end+n_ic) = csense_ic;
0198 
0199 % vobj = v_ic*b_ic + Sj ub(j)*uu(j)
0200 A_bl = [A_bl; c' sparse(1,n_int+m) -ub' sparse(1,n_int) -b_ic' sparse(1,2*n_m)];
0201 b_bl = [b_bl; 0];
0202 csense_bl(end+1) = 'E';
0203 
0204 % v(j) <= ub(j)*y(j) j in selected
0205 A_bl = [A_bl; Iint sparse(diag(-ub(sel_int == 1))) sparse(n_int,m+n+n_int+n_ic+2*n_m)];
0206 b_bl = [b_bl; zeros(n_int,1)];
0207 csense_bl(end+1:end+n_int) = 'L';
0208 
0209 % v(j) <= ub(j) j not in selected or inner constraints
0210 A_bl = [A_bl; Inint sparse(n_nint,m+n+2*n_int+n_ic+2*n_m)];
0211 b_bl = [b_bl; ub(~sel_int & ~sel_ic)];
0212 csense_bl(end+1:end+n_nint) = 'L';
0213 
0214 % Dual for inner constraints
0215 A_bl = [A_bl; sparse(n_ic,n+n_int) S(:,ind_ic)' Iic sparse(n_ic,n_int) speye(n_ic) sparse(n_ic,2*n_m)];
0216 b_bl = [b_bl; c(ind_ic)];
0217 csense_bl(end+1:end+n_ic) = 'G';
0218 
0219 % Dual bound for n_int ub's
0220 A_bl = [A_bl; sparse(n_nint,n+n_int) S(:,sel_int == 0 & sel_ic == 0)' Inint sparse(n_nint,n_int+n_ic+2*n_m)];
0221 b_bl = [b_bl; zeros(n_nint,1)];
0222 csense_bl(end+1:end+n_nint) = 'G';
0223 
0224 % Dual bound for selected ub's
0225 A_bl = [A_bl; sparse(n_int,n) -H*speye(n_int) S(:,sel_int == 1)' Iint sparse(n_int,n_int+n_ic+2*n_m)];
0226 b_bl = [b_bl; -H*ones(n_int,1)];
0227 csense_bl(end+1:end+n_int) = 'G';
0228 
0229 % Dual bound for selected <= 0's
0230 A_bl = [A_bl; sparse(n_int,n) H*speye(n_int) S(:,sel_int == 1)' sparse(n_int,n) speye(n_int) sparse(n_int,n_ic+2*n_m)];
0231 b_bl = [b_bl; zeros(n_int,1)];
0232 csense_bl(end+1:end+n_int) = 'G';
0233 
0234 % Dual bound for selected >= 0's
0235 A_bl = [A_bl; sparse(n_int,n) -H*speye(n_int) S(:,sel_int == 1)' sparse(n_int,n) speye(n_int) sparse(n_int,n_ic+2*n_m)];
0236 b_bl = [b_bl; zeros(n_int,1)];
0237 csense_bl(end+1:end+n_int) = 'L';
0238 
0239 % Create matchings for reversible rxns
0240 rev_match = sparse(n_intr,n_int);
0241 sel_fw_rxn = ones(1,n_int);
0242 for i = 1:n_intr
0243     rev_match(i,rev_int(i,1)) = 1;
0244     rev_match(i,rev_int(i,2)) = -1;
0245     sel_fw_rxn(rev_int(i,1)) = 0;
0246 end
0247 A_bl = [A_bl; sparse(n_intr,n) rev_match sparse(n_intr,m+n+n_int+n_ic+2*n_m)];
0248 b_bl = [b_bl; zeros(n_intr,1)];
0249 csense_bl(end+1:end+n_intr) = 'E';
0250 
0251 % Limit maximum number of deletions
0252 A_bl = [A_bl; sparse(1,n) -sel_fw_rxn sparse(1,m+n+n_int+n_ic+2*n_m)];
0253 b_bl = [b_bl; -sum(sel_fw_rxn) + K];
0254 csense_bl(end+1) = ksense;
0255 
0256 % Slack constraint for measured values
0257 if (n_m > 0)
0258     tmp1 = [Im; Im];
0259     tmp2 = [speye(n_m) sparse(n_m,n_m); sparse(n_m,n_m) -speye(n_m)];
0260     A_bl = [A_bl; tmp1 sparse(2*n_m,m+n+2*n_int+n_ic) tmp2];
0261     b_bl = [b_bl; b_m; b_m];
0262     csense_bl(end+1:end+n_m) = 'G';
0263     csense_bl(end+1:end+n_m) = 'L';
0264 end
0265 
0266 % Add constraint to select a new deletion set (with at least one new deletion)
0267 if (~isempty(sel_int_prev))
0268     [tmp,n_exko] = size(sel_int_prev);
0269     for i = 1:n_exko
0270         A_bl = [A_bl; sparse(1,n) sparse(sel_int_prev(sel_int==1,i)') sparse(1,n+n_int+m+n_ic+2*n_m)];
0271         b_bl(end+1) = 1;
0272         csense_bl(end+1) = 'G';
0273     end
0274 end
0275 
0276 % Construct problem structure
0277 bilevelMILPproblem.A = A_bl;
0278 bilevelMILPproblem.b = b_bl;
0279 bilevelMILPproblem.c = c_bl;
0280 bilevelMILPproblem.csense = csense_bl;
0281 bilevelMILPproblem.lb = lb_bl;
0282 bilevelMILPproblem.ub = ub_bl;
0283 bilevelMILPproblem.vartype = vartype_bl;
0284 bilevelMILPproblem.contSolInd = sel_cont_sol;
0285 bilevelMILPproblem.intSolInd = sel_int_sol;

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