Case study Stochastic Utility Problem (see Formal Problem Statement) in MATLAB Environment is solved with mpsg_solver PSG function.
Main MATLAB code is in file CS_Stochastic_Utility_Problem_Scen_Lin_Blocks_Toolbox.m.
Data are saved in file CS_Utility_data.zip:
Let us describe the main operations. To run case study you need to do the following main steps:
In file CS_Stochastic_Utility_Problem_Scen_Lin_Blocks_Toolbox.m:
Reading and checking problem data:
fid = fopen('l500.ins','r');
if fid<0, disp(' '); disp(' DATA_file does not exist'); break;
end
Read n:
rline=fgetl(fid);
if rline == -1, disp(' '); disp(['Error in ', file_name, ': Keyword "n=" is not found']); terminate %End %quit % Stop
else; rline=strtrim(lower(rline)); rl=length(rline);
if(rl>0); nw = sscanf(rline, '%i'); n=nw(1);
end
end
Read V_k, S_k:
k=0;
while true, rline=fgetl(fid);
if rline == -1, break;
else; rl=length(rline);
if(rl>2), [v(k+1),s(k+1)] = strread(rline); k=k+1; else, break;
end
end
end
fclose(fid);
Read scenarios for scenarios problem statement:
ps = fopen('Utility_senario_forL500_Short.txt', 'r');
J=0; tic;
while true, rline=fgetl(ps);
if rline == -1, break; else; J=J+1; end
end
frewind(ps); disp(['counting scenarios ',num2str(toc)]);
matrix_scen(J,n)=0; l=0; tic;
while true, rline=fgetl(ps);
if rline == -1, break;
else; rl=length(rline); l=l+1;
ip = findstr(rline(8:14),':');
matrix_scen(l,1:n)=strread(rline(8+ip:rl), '%f' )';
end
end
disp(['reading scenarios ',num2str(toc)]);
Make pmatrix_scenario:
pmatrix_scen((n+1)*J,3)=0.; id=0; tic;
for j=1:J
for i=1:n; id=id+1;
pmatrix_scen(id,1) = j;
pmatrix_scen(id,2) = i;
pmatrix_scen(id,3) = matrix_scen(j,i);
end
id=id+1;
pmatrix_scen(id,1) = j;
pmatrix_scen(id,2) = n+j; % variable T_j
pmatrix_scen(id,3) = -1;
end
Spmatrix_scen = sparse(pmatrix_scen(:,1),pmatrix_scen(:,2),pmatrix_scen(:,3));
clear matrix_scen pmatrix_scen;
disp(['filling pmatrix of scenarios ',num2str(toc)]);
Pack matrices and point:
variables_X='x1'; variables_T='T_1';
for j=2:n; variables_X =[variables_X, char(9),'x',int2str(j)]; end
for j=2:J; variables_T =[variables_T, char(9),'T_',int2str(j)]; end
count=1;
toolboxstruc_arr(count) = tbpsg_pmatrix_pack('pmatrix_scen', Spmatrix_scen, [variables_X char(9) variables_T]);
matrix_block_k(1,2)=0.;
for l=1:k
count=count+1;
toolboxstruc_arr(count) = tbpsg_matrix_pack(['matrix_',int2str(l)], -s(l),['T_1'],v(l));
end
avg_matrices=['block_',int2str(J),'(matrix_1)'];
for l=2:k; avg_matrices=[avg_matrices ', block_',int2str(J),'(matrix_',int2str(l),')'];
end
count=count+1;
toolboxstruc_arr(count) = tbpsg_matrix_pack('matrix_budget', ones(1,n),variables_X);
count=count+1;
toolboxstruc_arr(count) = tbpsg_point_pack('point_lowerbounds',zeros(1,n),variables_X);
Form problem statement:
problem_statement=sprintf('%s\n',...
'problem: problem_Utility_Scen_Block, type = minimize',...
'objective: objective_avg_multi, linearize = 0',...
[' avg_max_risk_(', avg_matrices ,')'],...
'constraint: constraint_linear_mult, lower_bound = 0, upper_bound = 0',...
' linearmulti_1(pmatrix_scen)',...
'constraint: constraint_budget, lower_bound = 1, upper_bound = 1',...
' linear_budget(matrix_budget)',...
'box_of_variables: lowerbounds = 0',...
'solver: car, precision = 7',...
'');
tic
%Uncomment the following line to open the problem in Toolbox Window:
%tbpsg_toolbox(problem_statement,toolboxstruc_arr);
Optimize problem:
[solution_str, outargstruc_arr] = tbpsg_run(problem_statement, toolboxstruc_arr); |
Display solution and optimal point:
disp(solution_str);
optimal_point_data = tbpsg_optimal_point_data(solution_str, outargstruc_arr);
optimal_point_vars = tbpsg_optimal_point_vars(solution_str, outargstruc_arr);
Problem: solution_status = optimal
Timing: Data_loading_time = 0.02, Preprocessing_time = 0.02, Solving_time = 0.02
Variables: optimal_point = point_problem_utility_scen_block
Objective: objective_avg_multi = -8.707287523413
Constraint: constraint_linear_mult = vector_constraint_linear_mult, [6.661338147751e-016]
Constraint: constraint_budget = 1.000000000000e+000, [-2.220446049250e-016]
Function: avg_max_risk_(block_13(matrix_1), block_13(matrix_2), block_13(matrix_3), block_13(matrix_4), block_13(matrix_5), block_13(matrix_6), block_13(matrix_7), block_13(matrix_8), block_13(matrix_9), block_13(matrix_10), block_13(matrix_11)) = -8.707287523413e+000
Function: linearmulti_1(pmatrix_scen) = vector_constraint_linear_mult
Function: linear_budget(matrix_budget) = 1.000000000000e+000