0001 function [statusOK, invalidConstraints, invalidVars, objective] = verifyCobraProblem(XPproblem, x, tol, verbose)
0002
0003
0004
0005
0006
0007
0008
0009
0010
0011
0012
0013
0014
0015
0016
0017
0018
0019
0020
0021
0022
0023
0024
0025
0026
0027
0028
0029
0030
0031
0032
0033
0034
0035
0036
0037
0038 if nargin < 3
0039 tol = 1e-8;
0040 end
0041 if nargin < 4
0042 verbose = true;
0043 end
0044
0045
0046 validQP = false;
0047 validMI = false;
0048 objective = [];
0049 statusOK = 1;
0050
0051
0052 if ~isfield(XPproblem, 'A')
0053 display('Required field A not found');
0054 statusOK = -1;
0055 return;
0056 elseif any(isnan(XPproblem.A))
0057 [r c] = find(isnan(XPproblem.A));
0058 strCoords = '';
0059 for i=1:length(r)
0060 strCoords = [strCoords ' ' num2str(r(i)) ',' num2str(c(i))];
0061 end
0062 display(['NaN present in A matrix at' strCoords '.']);
0063 statusOK = -1;
0064 return;
0065 end
0066 [nconstraints, nvars] = size(XPproblem.A);
0067
0068
0069
0070 if ~isfield(XPproblem, 'b')
0071 display('Required field b not found');
0072 statusOK = -1;
0073 return;
0074 elseif any(isnan(XPproblem.b))
0075 r= find(isnan(XPproblem.b));
0076 strCoords = '';
0077 for i=1:length(r)
0078 strCoords = [strCoords ' ' num2str(r(i)) ','];
0079 end
0080 display(['NaN present in b vector at' strCoords '.']);
0081 statusOK = -1;
0082 return;
0083 end
0084 if any(size(XPproblem.b) ~= [nconstraints, 1])
0085 display('Wrong size b vector');
0086 statusOK = -1;
0087 return;
0088 end
0089
0090
0091 if ~isfield(XPproblem, 'csense')
0092 display('Required field csense not found');
0093 statusOK = -1;
0094 return;
0095 end
0096 if length(XPproblem.csense) ~= nconstraints
0097 display('Wrong size csense vector');
0098 statusOK = -1;
0099 return;
0100 end
0101 if size(XPproblem.csense,2) ~= 1
0102 display('Csense should be a column vector')
0103 statusOK = -1;
0104 return;
0105 end
0106 invalidCsense = find(~ismember(XPproblem.csense,['E','L','G']));
0107 if invalidCsense
0108 fprintf('Invalid csense entry(s) at %s\n', num2str(invalidCsense));
0109 statusOK = -1;
0110 return;
0111 end
0112
0113
0114 if ~isfield(XPproblem,'lb')
0115 diplay('Required field lb not found');
0116 statusOK = -1;
0117 return;
0118 elseif any(isnan(XPproblem.lb))
0119 r= find(isnan(XPproblem.lb));
0120 strCoords = '';
0121 for i=1:length(r)
0122 strCoords = [strCoords ' ' num2str(r(i)) ','];
0123 end
0124 display(['NaN present in lb vector at' strCoords '.']);
0125 statusOK = -1;
0126 return;
0127 end
0128 if any(size(XPproblem.lb) ~= [nvars, 1])
0129 display('Wrong size lb vector');
0130 statusOK = -1;
0131 return;
0132 end
0133
0134
0135 if ~isfield(XPproblem,'ub')
0136 diplay('Required field ub not found');
0137 statusOK = -1;
0138 return;
0139 elseif any(isnan(XPproblem.ub))
0140 r= find(isnan(XPproblem.ub));
0141 strCoords = '';
0142 for i=1:length(r)
0143 strCoords = [strCoords ' ' num2str(r(i)) ','];
0144 end
0145 display(['NaN present in ub vector at' strCoords '.']);
0146 statusOK = -1;
0147 return;
0148 end
0149 if any(size(XPproblem.ub) ~= [nvars, 1])
0150 display('Wrong size ub vector');
0151 statusOK = -1;
0152 return;
0153 end
0154
0155 if any(XPproblem.ub<XPproblem.lb)
0156 fprintf('Invalid lb/ub at %s\n', num2str(find(XPproblem.ub<XPproblem.lb)));
0157 statusOK = -1;
0158 return;
0159 end
0160
0161
0162 if ~isfield(XPproblem,'c')
0163 diplay('Required field c not found');
0164 statusOK = -1;
0165 return;
0166 elseif any(isnan(XPproblem.c))
0167 r= find(isnan(XPproblem.c));
0168 strCoords = '';
0169 for i=1:length(r)
0170 strCoords = [strCoords ' ' num2str(r(i)) ','];
0171 end
0172 display(['NaN present in c vector at' strCoords '.']);
0173 statusOK = -1;
0174 return;
0175 end
0176 if any(size(XPproblem.c) ~= [nvars, 1])
0177 display('Wrong size c vector');
0178 statusOK = -1;
0179 return;
0180 end
0181
0182 validLP = true;
0183
0184 if isfield(XPproblem,'F')
0185 [nRows nCols] = size(XPproblem.F);
0186 if any(isnan(XPproblem.F))
0187 [r c]= find(isnan(XPproblem.F));
0188 strCoords = '';
0189 for i=1:length(r)
0190 strCoords = [strCoords ' ' num2str(r(i)) ',' num2str(c(i))];
0191 end
0192 display(['NaN present in F matrix at' strCoords '.']);
0193 statusOK = -1;
0194 return;
0195 end
0196 if nRows ~= nCols
0197 display('F matrix not square');
0198 statusOK = -1;
0199 elseif nRows ~= nvars
0200 display('Wrong size F matrix');
0201 statusOK = -1;
0202 else
0203 validQP = true;
0204 end
0205 end
0206
0207 if isfield(XPproblem,'vartype')
0208 if all(size(XPproblem.vartype) == [nvars,1])
0209 invalidVartype = find(~ismember(XPproblem.vartype,['C','I','B']));
0210 if isempty(invalidVartype)
0211 validMI = true;
0212 else
0213 fprintf('Invalid vartype entry(s) at %s\n', num2str(invalidVartype));
0214 statusOK = -1;
0215 end
0216 else
0217 display('Wrong size vartype vector');
0218 statusOK = -1;
0219 end
0220 vartype = XPproblem.vartype;
0221 if any(floor(XPproblem.ub(vartype == 'I' | vartype == 'B') + tol) < ceil(XPproblem.lb(vartype =='I' | vartype == 'B') - tol))
0222 display('Integer or binary variables lb to ub range does not contain an integer');
0223 validMI = false;
0224 statusOK = -1;
0225 end
0226 if any(XPproblem.lb(vartype == 'B') ~= 0)
0227 display('Binary variables have lower bound not equal to zero. This is inconsistent');
0228 validMI = false;
0229 statusOK = -1;
0230 end
0231 if any(XPproblem.ub(vartype == 'B') ~= 1)
0232 display('Binary variables have upper bound not equal to one. This is inconsistent');
0233 validMI=false;
0234 statusOK = -1;
0235 end
0236 end
0237
0238 if verbose
0239 if validLP
0240 display('Valid LP problem');
0241 else
0242 display('Invalid LP problem');
0243 end
0244 if validMI && validLP
0245 display('Valid MILP problem');
0246 else
0247 display('Invalid MILP problem');
0248 end
0249 if validQP
0250 display('Valid QP problem');
0251 else
0252 display('Invalid QP problem');
0253 end
0254 if validMI && validQP
0255 display('Valid MIQP problem');
0256 else
0257 display('Invalid MIQP problem');
0258 end
0259 if ~validLP&&~validQP
0260 return;
0261 end
0262 end
0263
0264
0265 if nargin >= 2 && ~isempty(x)
0266 validX = true;
0267 validXMI = false;
0268 if any(size(x)~=[nvars,1])
0269 display('Wrong size x vector');
0270 statusOK = 0;
0271 return;
0272 end
0273 if any(isnan(x))
0274 r= find(isnan(x));
0275 strCoords = '';
0276 for i=1:length(r)
0277 strCoords = [strCoords ' ' num2str(r(i)) ','];
0278 end
0279 display(['NaN present in x vector at' strCoords '.']);
0280 statusOK = -1;
0281 return;
0282 end
0283 invalidConstraints = zeros(nconstraints,1);
0284 invalidVars = zeros(nvars,1);
0285 if any(x > XPproblem.ub + tol)
0286 invalidVars(x > XPproblem.ub + tol) = 1;
0287 display('Upper bound violation')
0288 statusOK = 0;
0289 end
0290 if any(x < XPproblem.lb - tol)
0291 invalidVars(x < XPproblem.lb - tol) = 1;
0292 display('Lower bound violation')
0293 statusOK = 0;
0294 end
0295 product = XPproblem.A*x;
0296
0297 if any(abs(product(XPproblem.csense == 'E') - XPproblem.b(XPproblem.csense == 'E')) > tol)
0298 invalidConstraints(abs(product(XPproblem.csense == 'E') - XPproblem.b(XPproblem.csense == 'E')) > tol) = 1;
0299 display('Equality constraint off');
0300 validX = false;
0301 statusOK = 0;
0302 end
0303 if any(product(XPproblem.csense == 'L') > XPproblem.b(XPproblem.csense == 'L') + tol)
0304 invalidConstraints(product(XPproblem.csense == 'L') > XPproblem.b(XPproblem.csense == 'L') + tol) = 1;
0305 display('L constraint off');
0306 validX = false;
0307 statusOK = 0;
0308 end
0309 if any(product(XPproblem.csense == 'G') < XPproblem.b(XPproblem.csense == 'G') - tol)
0310 invalidConstraints(product(XPproblem.csense == 'G') < XPproblem.b(XPproblem.csense == 'G') - tol) = 1;
0311 display('G constraint off');
0312 validX = false;
0313 statusOK = 0;
0314 end
0315
0316
0317 if isfield(XPproblem, 'vartype')
0318 validXMI = true;
0319 if(abs( x(vartype == 'I' | vartype == 'B') - round(x(vartype == 'I' | vartype == 'B'))) > tol)
0320 display('Integer constraint off')
0321 validXMI = false;
0322 statusOK = 0;
0323 end
0324 end
0325 if validX
0326 if validXMI
0327 display('Valid x vector for MIXP problem');
0328 else
0329 display('Valid x vector for XP problem');
0330 end
0331 end
0332
0333 if validQP
0334 objective = (1/2)*x'*XPproblem.F*x + XPproblem.c'*x;
0335 elseif validLP
0336 objective = XPproblem.c'*x;
0337 end
0338 invalidConstraints = find(invalidConstraints);
0339 invalidVars = find(invalidVars);
0340 end
0341
0342
0343
0344
0345