%This case study considers the problem of optimal selection of tests %subject to several constraints on available resources %(e.g. money, times, and people). There are no partial tests: each test %is assumed to be either conducted or not conducted. If each resource %estimate is assumed to be accurate, then the problem of optimal selection %of tests is formulated as Deterministic Linear Programming Assignment %model with boolean decision variables. To take into account uncertainty %in resource estimates two models are used: robust model and the %stochastic model. The Robust model conservatively increases the need %in each resource by 20% of its average consumption by 20% %largest consumers. The Stochastic model is based on the assumption that %resource consumption by each test is independent normally %distributed random value. clc clear; global CDF_table_0_1; global PDF_table_0_1; %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'); RHS=RHS(:,2); vector_RHS=RHS; %pack vector with right-hand side bounds toolboxstruc_arr(count) = tbpsg_vector_pack('vector_right', vector_RHS); count = count+1; %load tables of Normal CDF and PDF load('NormalDistribution'); %----------------- %Solve Problem 1 | %----------------- %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); %----------------- %Solve Problem 2 %----------------- %generate and pack matrix_c1, ..., matrix_c20 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 clear problem_statement; 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'); clear outargstruc_arr ; %call solver [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); %-------------------------------- %Solve Problem 3 with tbpsg_run %-------------------------------- %Build 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 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 matrix 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 clear problem_statement; 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'); clear outargstruc_arr %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| %----------------------------------------- probab_upper = 0.2; % upper bound on probability X = X_DET; %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); %gradient=gradient+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'); % 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 %Calculate probability for improved solution prob_opt = fnext_prmulti_pen_g_ni_m(point_variables, y_opt, mu_TETA,var_TETA,size(mu_TETA,1), 0); %output solutions of Problem 1, Problem 2, Problem 3 (PSG solver), and %Problem 3 (Heuristic Algorithm) fprintf('\nSolution of Problem 1:\n %s\n', solution_str1); fprintf('\nSolution of Problem 2:\n %s\n', solution_str2); fprintf('\nSolution of Problem Problem 3 (PSG solver):\n %s\n', solution_str3); fprintf('\nProblem 3 (Heuristic Algorithm): oblective = %12.5f\tprobability = %12.5f\n\n', opt_obj, prob_opt); fprintf('\nOptimal Points for Problem 1, Problem 2, Problem 3 (PSG solver), and Problem 3 (Heuristic Algorithm)'); fprintf('\n\n Component\tProblem 1\tProblem 2\tProblem 3 (PSG solver)\tProblem 3 (Heuristic Algorithm)\n'); for j=1:1:length(point_variables) fprintf('%5s\t%9.0f\t%9.0f\t%15.0f\t%20.0f\n', point_variables{j},X_DET(j),X_2(j),X_3(j), y_opt(j)); end %=======================================================================| %American Optimal Decisions, Inc. Copyright | %Copyright ©American Optimal Decisions, Inc. 2007-2016. | %American Optimal Decisions (AOD) retains copyrights to this material. | % | %Permission to reproduce this document and to prepare derivative works | %from this document for internal use is granted, provided the copyright | %and “No Warranty” statements are included with all reproductions | %and derivative works. | % | %For information regarding external or commercial use of copyrighted | %materials owned by AOD, contact AOD at support@aorda.com. | %=======================================================================|