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
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;