Solution in MATLAB Toolbox
Contents
Problem 1. Deterministic Linear Programming Assignment Model
Case study Optimal Tests Selection (see Formal Problem Statement) in MATLAB Environment is solved with tbpsg_run PSG function.
Main MATLAB code is in file CS_Optimal_Tests_Selection_Toolbox.m.
Data are saved in file Optimal_Tests_Selection_data.mat:
matrix_obj |
Values of tests for objective functions (CS.1), (CS.4), (CS.7). |
matrix_multilin |
Coefficients of resources consumption for constraints (CS.2), (CS.5) . These coefficients are used for building Probability Exceeding Penalty for Gain Multiple Normal Independent function in constraint (CS.8). |
vector_RHS |
Upper bounds for constraints (CS.2), (CS.5) |
types_arr |
Types of variables for constraints for the budget constraint (CS.3), (CS.6), (CS.9). |
Let us describe the main operations. To run case study you need to do the following main steps:
In file CS_Optimal_Tests_Selection.m.
Prepare data in PSG format for solving problems in PSG General (text) format in MATLAB:
Set parameters values:
C_sigma = 0.2;% is used for calculation of st.dev.: st.dev. = C_sigma*average
w = 0; %threshold
Load data:
load('.\Optimal_Tests_Selection_data.mat');
Pack matrix for building objective:
count = 1;
clear toolboxstruc_arr;
toolboxstruc_arr(count) = tbpsg_matrix_pack('matrix_obj', matrix_obj, header_matrix_obj);
count = count+1;
Pack point with types:
toolboxstruc_arr(count) = tbpsg_point_pack('point_types', types_arr, header_point);
count = count+1;
Pack matrix for multilinear constraints:
toolboxstruc_arr(count) = tbpsg_matrix_pack('matrix_multilin', matrix_multilin, header_matrix_multilin);
count = count+1;
Create array of columns names in matrix_multilin:
matrices_names = strread(header_matrix_multilin, '%s');
Pack vector with right-hand side bounds:
RHS=RHS(:,2);
vector_RHS=RHS;
toolboxstruc_arr(count) = tbpsg_vector_pack('vector_right', vector_RHS);
count = count+1;
Load tables of Normal CDF and PDF:
load('NormalDistribution');
Problem 1. Deterministic Linear Programming Assignment Model (CS.1) - (CS.3)
Generate statement for Problem 1:
problem_statement = sprintf('%s\n',...
'maximize',...
'linear(matrix_obj)',...
'Constraint: <= vector_right',...
'linearmulti(matrix_multilin)',...
'Box: lowerbounds = 0, upperbounds = 1, types = 1',...
'Solver: init_point = 0.9');
Optimize problem:
[solution_str1, outargstruc_arr] = tbpsg_run(problem_statement, toolboxstruc_arr); |
Extract optimal point for Problem 1:
point_variables = tbpsg_optimal_point_vars(solution_str1, outargstruc_arr);
X_DET = tbpsg_optimal_point_data(solution_str1, outargstruc_arr);
Problem 2. Robust Model (CS.4) - (CS.6)
For solving the Problem 2 divide 20-row matrix for multilinear function into 20 one-row matrices for linear functions:
for i = 1:1:20
clear matrix_c matrix_name;
matrix_c = matrix_multilin(i,:);
matrix_c(1) = 1;
matrix_name = ['matrix_c' int2str(i)];
options.create = 1;
toolboxstruc_arr(count) = tbpsg_matrix_pack(matrix_name, matrix_c, header_matrix_multilin,[],[],options);
count = count+1;
end
toolboxstruc_arr(count) = tbpsg_vector_pack('vector_upper_bound',[800,600,6,5,6,5,6,5,6,5,6,6,9,5,20,9,20,9,20,200]);
Generate statement for Problem 2:
problem_statement = sprintf('%s\n',...
'maximize',...
'linear(matrix_obj)',...
'MultiConstraint: <= vector_upper_bound',...
'linear(matrix_c1,...,matrix_c20)',...
'+0.2*cvar_comp_pos(0.8, matrix_c1,...,matrix_c20)',...
'Box: >= 0, <= 1, types = 1',...
'Solver: CAR, precision = 9');
Optimize problem:
[solution_str2, outargstruc_arr] = tbpsg_run(problem_statement, toolboxstruc_arr); |
Extract optimal point for Problem 2:
point_variables2 = tbpsg_optimal_point_vars(solution_str2, outargstruc_arr);
X_2 = tbpsg_optimal_point_data(solution_str2, outargstruc_arr);
Problem 3. Stochastic Model (CS.7) - (CS.9)
Build mean and variance matrices needed for modeling normal distribution.
Build mean matrix "mu_TETA":
mu_TETA = [];
for j = 1:1:size(matrix_multilin,2)
disp(j)
if strcmpi(matrices_names{j},'scenario_benchmark')
mu_TETA = [ mu_TETA RHS + matrix_multilin(:,j)];
else
mu_TETA = [ mu_TETA matrix_multilin(:,j)];
end
end
Build variance matrix "var_TETA":
A = (C_sigma*mu_TETA(:,2:size(matrix_multilin,2)-2)).^2;
var_TETA = [mu_TETA(:,1) A mu_TETA(:,size(matrix_multilin,2)-1) zeros(size(matrix_multilin,1),1)];
Pack matrices mu_TETA and var_TETA:
toolboxstruc_arr(count) = tbpsg_matrix_pack('matrix_mu_TETA', mu_TETA, header_matrix_multilin);
count = count+1;
toolboxstruc_arr(count) = tbpsg_matrix_pack('matrix_var_TETA', var_TETA, header_matrix_multilin);
count = count+1;
Generate statement for Problem 3:
problem_statement = sprintf('%s\n',...
'maximize',...
'linear(matrix_obj)',...
'Constraint: <= vector_right',...
'linearmulti(matrix_multilin)',...
'Constraint:<= 0.2',...
'prmulti_pen_ni_g(0.0, matrix_mu_TETA, matrix_var_TETA)',...
'Box: >= 0, <= 1, types = 1',...
'Solver: init_point = 0.9');
Optimize problem:
[solution_str3, outargstruc_arr] = tbpsg_run(problem_statement, toolboxstruc_arr); |
Extract optimal point for Problem 3:
point_variables3 = tbpsg_optimal_point_vars(solution_str3, outargstruc_arr);
X_3 = tbpsg_optimal_point_data(solution_str3, outargstruc_arr);
Solve Problem 3 with Heuristic Algorithm.
Rebuild matrix mu_TETA:
mu_TETA = [];
headerpoint{1} ='scenario_benchmark';
for j = 2:1:size(matrix_multilin,2)
if strcmpi(matrices_names{j},'scenario_benchmark')
mu_TETA = [RHS + matrix_multilin(:,j) mu_TETA ];
else
if ~strcmpi(matrices_names{j},'scenario_probability')
headerpoint{j}=matrices_names(j);
mu_TETA = [ mu_TETA matrix_multilin(:,j)];
end
end
end
Rebuild matrix var_TETA:
A = (C_sigma*mu_TETA(:,2:size(mu_TETA,2))).^2;
var_TETA = [0.001*ones(size(mu_TETA,1),1) A ];
Calculate objective value and probability for deterministic solution:
obj = matrix_obj(2:size(matrix_multilin,2)-2);
prob = fnext_prmulti_pen_g_ni_m(point_variables, X, mu_TETA,var_TETA,size(mu_TETA,1), 0);
prob_DET = prob;
obj1 = 0;
for i = 1:1:size(X,2)
obj1 = obj1 + obj(i)*X(i);
end
obj_DET = obj1;
Order deterministic solution with respect to values of X and values of ratio obective/gradient:
X_1 =[]; % set of components that equal 1
X_0 = []; % set of components that equal 0
for i = 1:1:size(X,2)
if X(1,i) == 1
X_1 =[X_1 i];
else
X_0 =[X_0 i];
end
end
y = X;
X_DET = X;
gradient = grext_prmulti_pen_g_ni_m(point_variables, X, mu_TETA,var_TETA,size(mu_TETA,1), 0);
for i = 1:1:size(X_1,2)
obj_per_grad_1(i) = obj(X_1(i))/gradient(X_1(i));
end
for i = 1:1:size(X_0,2)
obj_per_grad_0(i) = obj(X_0(i))/gradient(X_0(i));
end
[B_1,IX_1] = sort(obj_per_grad_1, 'descend');
[B_0,IX_0] = sort(obj_per_grad_0, 'descend');
Transfer deterministic solution to solution that satisfies to probabilistic constraint (feasible solution):
j = size(IX_1,2);
while probab_upper < prob
y(X_1(IX_1(j))) = 0;
prob = fnext_prmulti_pen_g_ni_m(point_variables, y, mu_TETA,var_TETA,size(mu_TETA,1), 0);
j = j-1;
end
X = y;
Order feasible solution with respect to values of X and values of ratio obective/gradient.
while true
X_1 =[];
X_0 = [];
for i = 1:1:size(X,2)
if X(1,i) == 1
X_1 =[X_1 i];
else
X_0 =[X_0 i];
end
end
y_opt = X;
opt_obj = 0;
for i = 1:1:size(X,2)
opt_obj = opt_obj + obj(i)*X(i);
end
gradient = grext_prmulti_pen_g_ni_m(point_variables, X, mu_TETA,var_TETA,size(mu_TETA,1), 0);
clear obj_per_grad_0 obj_per_grad_1 IX_1 IX_0 B_0 B_1;
for i = 1:1:size(X_1,2)
obj_per_grad_1(i) = obj(X_1(i))/gradient(X_1(i));
end
for i = 1:1:size(X_0,2)
obj_per_grad_0(i) = obj(X_0(i))/gradient(X_0(i));
end
[B_1,IX_1] = sort(obj_per_grad_1, 'descend');
[B_0,IX_0] = sort(obj_per_grad_0, 'descend');
Search improvement for feasible solution:
j1 =1;
while true
y(X_0(IX_0(j1))) = 1;
prob = fnext_prmulti_pen_g_ni_m(point_variables, y, mu_TETA,var_TETA,size(mu_TETA,1), 0);
i = size(IX_1,2);
while true
if prob <= probab_upper
break;
end
y(X_1(IX_1(i))) = 0;
prob = fnext_prmulti_pen_g_ni_m(point_variables, y, mu_TETA,var_TETA,size(mu_TETA,1), 0);
i = i-1;
end
obj1 = 0;
for i = 1:1:size(y,2)
obj1 = obj1 + obj(i)*y(i);
end
if obj1 < opt_obj
break;
end
opt_obj = obj1;
y_opt = y;
j1 = j1+1;
end
if X == y_opt
break;
end
X = y_opt;
end
Problem: solution_status = optimal
Timing: Data_loading_time = 0.02, Preprocessing_time = 0, Solving_time = 0.11
Variables: optimal_point = point_problem_selection_multilp_assignment_model
Objective: objective_linear = 879
Constraint: constraint_multilin = vector_constraint_multilin, constraint_multilin_residual = -2.458781320000e-001
Function: linear_obj(matrix_obj) = 8.790000000000e+002
Function: linearmulti_matrix_multi(matrix_multilin) = vector_constraint_multilin
Problem: solution_status = optimal
Timing: Data_loading_time = 0.05, Preprocessing_time = 0.03, Solving_time = 0.3
Variables: optimal_point = point_problem_robust
Objective: objective_robust = 873
Constraint: constraint_c1 = 5.490571428571e+002, constraint_c1_residual = -2.509428571429e+002
Constraint: constraint_c2 = 5.973333333333e+002, constraint_c2_residual = -2.666666666667e+000
Constraint: constraint_c3 = 3.570135721429e+000, constraint_c3_residual = -2.429864278571e+000
Constraint: constraint_c4 = 3.131543452190e+000, constraint_c4_residual = -1.868456547810e+000
Constraint: constraint_c5 = 1.498959057619e+000, constraint_c5_residual = -4.501040942381e+000
Constraint: constraint_c6 = 0.000000000000e+000, constraint_c6_residual = -5.000000000000e+000
Constraint: constraint_lc7 = 1.097890644190e+000, constraint_lc7_residual = -4.902109355810e+000
Constraint: constraint_c8 = 1.047619047619e+000, constraint_c8_residual = -3.952380952381e+000
Constraint: constraint_c9 = 2.159539716190e-001, constraint_c9_residual = -5.784046028381e+000
Constraint: constraint_c10 = 1.047619047619e+000, constraint_c10_residual = -3.952380952381e+000
Constraint: constraint_c11 = 1.017723288476e+000, constraint_c11_residual = -4.982276711524e+000
Constraint: constraint_c12 = 0.000000000000e+000, constraint_c12_residual = -6.000000000000e+000
Constraint: constraint_c13 = 8.843824536010e+000, constraint_c13_residual = -1.561754639905e-001
Constraint: constraint_c14 = 5.767712269524e-001, constraint_c14_residual = -4.423228773048e+000
Constraint: constraint_c15 = 7.856743767581e+000, constraint_c15_residual = -1.214325623242e+001
Constraint: constraint_c16 = 5.767712269524e-001, constraint_c16_residual = -8.423228773048e+000
Constraint: constraint_c17 = 1.009667171823e+001, constraint_c17_residual = -9.903328281771e+000
Constraint: constraint_c18 = 1.927201308476e+000, constraint_c18_residual = -7.072798691524e+000
Constraint: constraint_c19 = 8.722849303562e+000, constraint_c19_residual = -1.127715069644e+001
Constraint: constraint_c20 = 1.322761904762e+002, constraint_c20_residual = -6.772380952381e+001
Function: linear_obj(matrix_obj) = 8.730000000000e+002
Function: linear_1(matrix_c1) = 5.310000000000e+002
Function: cvar_comp_pos_1(0.800000E+00,matrix_c1) = 9.028571428571e+001
Function: linear_2(matrix_c2) = 5.800000000000e+002
Function: cvar_comp_pos_2(0.800000E+00,matrix_c2) = 8.666666666667e+001
Function: linear_3(matrix_c3) = 3.407856825000e+000
Function: cvar_comp_pos_3(0.800000E+00,matrix_c3) = 8.113944821429e-001
Function: linear_4(matrix_c4) = 2.989200568000e+000
Function: cvar_comp_pos_4(0.800000E+00,matrix_c4) = 7.117144209524e-001
Function: linear_5(matrix_c5) = 1.430824555000e+000
Function: cvar_comp_pos_5(0.800000E+00,matrix_c5) = 3.406725130952e-001
Function: linear_6(matrix_c6) = 0.000000000000e+000
Function: cvar_comp_pos_6(0.800000E+00,matrix_c6) = 0.000000000000e+000
Function: linear_7(matrix_c7) = 1.047986524000e+000
Function: cvar_comp_pos_7(0.800000E+00,matrix_c7) = 2.495206009524e-001
Function: linear_8(matrix_c8) = 1.000000000000e+000
Function: cvar_comp_pos_8(0.800000E+00,matrix_c8) = 2.380952380952e-001
Function: linear_9(matrix_c9) = 2.061378820000e-001
Function: cvar_comp_pos_9(0.800000E+00,matrix_c9) = 4.908044809524e-002
Function: linear_10(matrix_c10) = 1.000000000000e+000
Function: cvar_comp_pos_10(0.800000E+00,matrix_c10) = 2.380952380952e-001
Function: linear_11(matrix_c11) = 9.714631390000e-001
Function: cvar_comp_pos_11(0.800000E+00,matrix_c11) = 2.313007473810e-001
Function: linear_12(matrix_c12) = 0.000000000000e+000
Function: cvar_comp_pos_12(0.800000E+00,matrix_c12) = 0.000000000000e+000
Function: linear_13(matrix_c13) = 8.573748297000e+000
Function: cvar_comp_pos_13(0.800000E+00,matrix_c13) = 1.350381195048e+000
Function: linear_14(matrix_c14) = 5.505543530000e-001
Function: cvar_comp_pos_14(0.800000E+00,matrix_c14) = 1.310843697619e-001
Function: linear_15(matrix_c15) = 7.589216958000e+000
Function: cvar_comp_pos_15(0.800000E+00,matrix_c15) = 1.337634047905e+000
Function: linear_16(matrix_c16) = 5.505543530000e-001
Function: cvar_comp_pos_16(0.800000E+00,matrix_c16) = 1.310843697619e-001
Function: linear_17(matrix_c17) = 9.740611634000e+000
Function: cvar_comp_pos_17(0.800000E+00,matrix_c17) = 1.780300421143e+000
Function: linear_18(matrix_c18) = 1.839601249000e+000
Function: cvar_comp_pos_18(0.800000E+00,matrix_c18) = 4.380002973810e-001
Function: linear_19(matrix_c19) = 8.463486388000e+000
Function: cvar_comp_pos_19(0.800000E+00,matrix_c19) = 1.296814577810e+000
Function: linear_20(matrix_c20) = 1.280000000000e+002
Function: cvar_comp_pos_20(0.800000E+00,matrix_c20) = 2.138095238095e+001
Problem: solution_status = optimal
Timing: Data_loading_time = 0.02, Preprocessing_time = 0, Solving_time = 0.28
Variables: optimal_point = point_problem_selection_multilp_assignment_model
Objective: objective_acc = 833
Constraint: constraint_multilin = vector_constraint_multilin, constraint_multilin_residual = -1.296978524000e+000
Constraint: constraint_prmulti_ni20 = 4.006307485245e-002, constraint_prmulti_ni20_residual = -1.599369251475e-001
Function: linear_obj(matrix_obj) = 8.330000000000e+002
Function: linearmulti_matrix_multi(matrix_multilin) = vector_constraint_multilin
Function: prmulti_pen_ni_g_20(0.000000E+00,matrix_mu_teta, matrix_var_teta) = 4.006307485245e-002
Optimal Points for Problem 1, Problem 2, Problem 3 (PSG solver), and Problem 3 (Heuristic Algorithm)
Component Problem 1 Problem 2 Problem 3 (PSG solver) Problem 3 (Heuristic Algorithm)
a 1 1 1 1
b 1 1 1 1
c 1 1 1 1
d 1 1 1 1
e 1 1 0 0
f 0 0 0 0
h 1 1 1 1
i 1 0 1 1
j 1 1 1 1
k 0 1 1 0
l 1 1 1 1
m 0 0 0 0
n 0 0 0 0
o 0 0 0 0
p 0 0 0 0
q 0 0 0 0
r 1 1 0 1
s 0 0 1 1
t 0 0 0 0
u 0 0 0 0
v 0 0 0 0