% This case study finds the best planting plan and insurance policy
% creating maximum profit under an acceptable risk exposure. 
% Farmers face uncertainties of crop yields and market prices affecting
% their profit. These uncertainties contribute to the profit risk.
% To hedge the risk, farmers can purchase insurance against the 
% uncertainties of yields and prices. There are several insurance policies
% on the market providing different protections with different prices. 
% Crop yields are affected by climate type and planting dates.
% Climate can be categorized into three phenomena: El Nino, La Nina and
% Neutral. Each phenomenon has different effects on the crop yields due to
% the different temperature and rainfall.
% In addition to crop yilds, the harvest price of crops is another
% important factor affecting farmers’ income. This factor is taken into
% consideration by using scenarios of the prices of crops, which were
% generated.
% The following three main types of crop insurance are considered: the
% Actual Production History crop % insurance (APH), 
% the Crop Revenue Coverage insurance (CRC), 
% and the Catastrophic Insurance Coverage (CAT). 
% APH assures a percentage of the farmers’ history yield. 
% CRC assures income by indemnifying farmers based on historical average
% yield and the market price. 
% CAT can be defined as an APH policy at 50% yield coverage with 55%
% price base election.
% We consider a model for planting and insurance by maximizing the total
% profit under a risk level measured by CVaR, VaR, 
% and Probability Exceeding Penalty.

clc
clear;

load Crop_Insurance_data_Toolbox.mat %toolboxstruc_arr
tbpsg_export_to_workspace(toolboxstruc_arr);
thistime = cputime;

% generate list of matrices names used in objective
str_ = ['matrix_obj_el ' 'matrix_obj_el ' 'matrix_obj_el ' 'matrix_obj_la ' 'matrix_obj_la ' 'matrix_obj_la ' 'matrix_obj_nt ' 'matrix_obj_nt ' 'matrix_obj_nt'];
matrix_obj = strread(str_,'%s');

% generate list of matrices names used in risk constraint
clear str_;
str_ = ['matrix_scenarios_el ' 'matrix_scenarios_el ' 'matrix_scenarios_el ' 'matrix_scenarios_la ' 'matrix_scenarios_la ' 'matrix_scenarios_la ' 'matrix_scenarios_nt ' 'matrix_scenarios_nt ' 'matrix_scenarios_nt'];
matrix_scen = strread(str_,'%s');

% generate list of fragments of risk functions with parameters used in problem
% description
clear str_;
str_ = ['cvar_risk_1(0.9, ' 'pr_pen_1(600, ' 'var_risk_1(0.9, ' 'cvar_risk_1(0.9, ' 'pr_pen_1(600, ' 'var_risk_1(0.9, ' 'cvar_risk_1(0.9, ' 'pr_pen_1(600, ' 'var_risk_1(0.9, '];
risk_finction =  strread(str_,'%s');

% generate lists of problem parameters
clear str_;
str_ = ['6000 ' '0.1 ' '925.226 ' '6000 ' '0.12 ' '737.146 ' '6000 ' '0.1 ' '-552.586 '];    
risk_upper = strread(str_, '%s');
clear str_;
str_ = ['Apr,16 ' 'Apr,23 ' 'May,1 ' 'May,8 ' 'May,15 ' 'May,22 ' 'May,29 ' 'Jun,5 ' 'Jun,12'];
planting_dates = strread(str_, '%s');
clear str_;
str_ = ['cotton ' 'peanut'];
crop_type = strread(str_, '%s');
clear str_;
str_ = ['No_insurance  ' 'CAT ' '65%APH ' '70%APH ' '75%APH ' '65%CRC  ' '70%CRC  ' '75%CRC  ' '80%CRC  ' '85%CRC  '  ];
insurance_policy  = strread(str_, '%s');


% generate the first problem statement
    str1 = sprintf('Problem: problem_crop_insurance, type = maximize'); 
    str2 = sprintf('Objective: objective_max_profit');
    str3 = sprintf('%s%s%s','linear_obj(',matrix_obj{1},')');
    str4 = sprintf('Constraint: constraint_acreage_cotton, lower_bound = 50, upper_bound = 50');
    str5 = sprintf('linear_acreage_cotton(matrix_acreage_cotton)');
    str6 = sprintf('Constraint: constraint_acreage_peanut, lower_bound = 50, upper_bound = 50');
    str7 = sprintf('linear_acreage_peanut(matrix_acreage_peanut)');
    str8 = sprintf('%s%s%s','Constraint: constraint_risk, upper_bound = ',risk_upper{1}, ', linearize = 1');
    str9 = sprintf('%s%s%s',risk_finction{1}, matrix_scen{1},')');
    str10 = sprintf('Constraint: constraint_acreage_xa, upper_bound = 0');
    str11 = sprintf('linearmulti_matrix_acreage_xa(matrix_acreage_xa)');
    str12 = sprintf('Constraint: constraint_cardinality_cotton, upper_bound = 1, linearize = 1');
    str13 = sprintf('cardn_pos_matrix_cardinality_cotton(0.000001, matrix_cardinality_cotton)');
    str14 = sprintf('Constraint: constraint_cardinlity_peanut, upper_bound = 1, linearize = 1');
    str15 = sprintf('cardn_pos_matrix_cardinality_peanut(0.000001, matrix_cardinality_peanut)');
    str16 = sprintf('Box_of_Variables: lowerbounds = point_lb, upperbounds = point_ub');
    str17 = sprintf('Solver: TANK, precision = 9, stages = 30');
    clear problem_statement;
    problem_statement = sprintf('%s\n', str1, str2, str3, str4, str5, str6, str7, str8, str9, str10, str11, str12, str13, str14, str15, str16, str17);

%Uncomment the following line to open the problem in Toolbox Window:
%tbpsg_toolbox(problem_statement,toolboxstruc_arr);    
    
% optimize the first problem      
    [solution_str, outargstruc_arr] = tbpsg_run(problem_statement, toolboxstruc_arr);

% extract optimal solution of the first problem 
point_variables = tbpsg_optimal_point_vars(solution_str, outargstruc_arr);
expected_profit(1) = tbpsg_objective(solution_str, outargstruc_arr); 
constr_data = tbpsg_constraints_data(solution_str, outargstruc_arr);
risk(1) = constr_data(3);

point_data = tbpsg_optimal_point_data(solution_str, outargstruc_arr);

for j=1:size(point_data,2)
   point_solution(j,1)=0;
   if abs(point_data(j)) > 0
      if log10(abs(point_data(j)))> -precision
         point_solution(j,1)= point_data(j); 
      end
   end
end


% begin loop over other problems    
    for i=2:1:9

% generate case-specific fragments of the next problem statement
        clear str3;
        str3 = sprintf('%s%s%s','linear_obj(',matrix_obj{i},')');
        clear str8;
        str8 = sprintf('%s%s%s','Constraint: constraint_risk, upper_bound = ',risk_upper{i}, ', linearize = 1');
        clear str9;
        str9 = sprintf('%s%s%s',risk_finction{i}, matrix_scen{i},')');
        clear str17;
        if i==2
           str17 = sprintf('Solver: VAN, precision = 6, stages = 30'); 
        end
        if i==3
           str17 = sprintf('Solver: VAN, precision = 6, stages = 30'); 
        end
         if i==4
           str17 = sprintf('Solver: TANK, precision = 9, stages = 30'); 
        end
        if i==5
           str17 = sprintf('Solver: VAN, precision = 6, stages = 20'); 
        end
        if i==6
           str17 = sprintf('Solver: VAN, precision = 9, stages = 30'); 
        end
        if i==7
           str17 = sprintf('Solver: TANK, precision = 6, stages = 20'); 
        end
        if i==8
           str17 = sprintf('Solver: VAN, precision = 6, stages = 20'); 
        end
        if i==9
           str17 = sprintf('Solver: VAN, precision = 6, stages = 30'); 
        end

% assemble and verify the current problem statement
        clear problem_statement;
        problem_statement = sprintf('%s\n', str1, str2, str3, str4, str5, str6, str7, str8, str9, str10, str11, str12, str13, str14, str15, str16, str17);
                
% optimize the current problem      
        [solution_str, outargstruc_arr] = tbpsg_run(problem_statement, toolboxstruc_arr);
        
% extract optimal solution of the current problem 
    point_variables = tbpsg_optimal_point_vars(solution_str, outargstruc_arr);
    expected_profit(i) = tbpsg_objective(solution_str, outargstruc_arr); 
    constr_data = tbpsg_constraints_data(solution_str, outargstruc_arr);
    risk(i) = constr_data(3);

    point_data = tbpsg_optimal_point_data(solution_str, outargstruc_arr);
    
    for j=1:size(point_data,2)
        point_solution(j,i)=0;
        if abs(point_data(j)) > 0
            if log10(abs(point_data(j)))> -precision
                point_solution(j,i)= point_data(j); 
            end
        end
    end

end
    
    fprintf('\nCPU time = %.2f%4s\n\n', cputime-thistime, 'sec');
    
% display and save tables with results of calculations
outfilename = 'CS_Crop_Insurance_Output.txt';
fid = fopen(outfilename, 'w');
fprintf('\n\n Table 1. Optimal Crop Production and Insurance Coverage for El Nino Climate Type');
fprintf(fid,'\n\n Table 1. Optimal Crop Production and Insurance Coverage for El Nino Climate Type');
fprintf('\n%78s', 'Planting area (in acres)'); 
fprintf(fid,'\n%78s', 'Planting area (in acres)'); 
fprintf('\n%10s%11s%15s%18s%6s%8s%6s', 'Variable', 'Crop Type', 'Planting Date', 'Insurance Policy', 'CVaR', 'PrPen', 'VaR');
fprintf(fid,'\n%10s%11s%15s%18s%6s%8s%6s', 'Variable', 'Crop Type', 'Planting Date', 'Insurance Policy', 'CVaR', 'PrPen', 'VaR');
for j=1:1:n_variables_principal
    if (point_solution(j,1)> 0) | (point_solution(j,2)> 0) | (point_solution(j,3)> 0)
        k1 = str2num(point_variables{j}(2))+1;
        k2 = str2num(point_variables{j}(3))+1;
        k3 = str2num(point_variables{j}(4))+1;
       fprintf('\n%9s%10s%14s%17s%10.2f%8.2f%8.2f', point_variables{j}, crop_type{k1}, planting_dates{k2}, insurance_policy{k3},  point_solution(j,1), point_solution(j,2), point_solution(j,3));  
       fprintf(fid,'\n%9s%10s%14s%17s%10.2f%8.2f%8.2f', point_variables{j}, crop_type{k1}, planting_dates{k2}, insurance_policy{k3},  point_solution(j,1), point_solution(j,2), point_solution(j,3));  
    end  
end

fprintf('\n\n Table 2. Optimal Crop Production and Insurance Coverage for La Nina Climate Type');
fprintf(fid,'\n\n Table 2. Optimal Crop Production and Insurance Coverage for La Nina Climate Type');
fprintf('\n%78s', 'Planting area (in acres)');
fprintf(fid,'\n%78s', 'Planting area (in acres)');
fprintf('\n%10s%11s%15s%18s%6s%8s%6s', 'Variable', 'Crop Type', 'Planting Date', 'Insurance Policy', 'CVaR', 'PrPen', 'VaR');
fprintf(fid,'\n%10s%11s%15s%18s%6s%8s%6s', 'Variable', 'Crop Type', 'Planting Date', 'Insurance Policy', 'CVaR', 'PrPen', 'VaR');
for j=1:1:n_variables_principal
    if (point_solution(j,4)> 0) | (point_solution(j,5)> 0) | (point_solution(j,6)> 0)
        k1 = str2num(point_variables{j}(2))+1;
        k2 = str2num(point_variables{j}(3))+1;
        k3 = str2num(point_variables{j}(4))+1;
       fprintf('\n%9s%10s%14s%17s%10.2f%8.2f%8.2f', point_variables{j}, crop_type{k1}, planting_dates{k2}, insurance_policy{k3},  point_solution(j,4), point_solution(j,5), point_solution(j,6));  
       fprintf(fid,'\n%9s%10s%14s%17s%10.2f%8.2f%8.2f', point_variables{j}, crop_type{k1}, planting_dates{k2}, insurance_policy{k3},  point_solution(j,4), point_solution(j,5), point_solution(j,6));  
    end  
end

fprintf('\n\n Table 3. Optimal Crop Production and Insurance Coverage for Neutral Climate Type');
fprintf(fid,'\n\n Table 3. Optimal Crop Production and Insurance Coverage for Neutral Climate Type');
fprintf('\n%78s', 'Planting area (in acres)'); 
fprintf(fid,'\n%78s', 'Planting area (in acres)');
fprintf('\n%10s%11s%15s%18s%6s%8s%6s', 'Variable', 'Crop Type', 'Planting Date', 'Insurance Policy', 'CVaR', 'PrPen', 'VaR');
fprintf(fid,'\n%10s%11s%15s%18s%6s%8s%6s', 'Variable', 'Crop Type', 'Planting Date', 'Insurance Policy', 'CVaR', 'PrPen', 'VaR');
for j=1:1:n_variables_principal
    if (point_solution(j,7)> 0) | (point_solution(j,8)> 0) | (point_solution(j,9)> 0)
        k1 = str2num(point_variables{j}(2))+1;
        k2 = str2num(point_variables{j}(3))+1;
        k3 = str2num(point_variables{j}(4))+1;
       fprintf('\n%9s%10s%14s%17s%10.2f%8.2f%8.2f', point_variables{j}, crop_type{k1}, planting_dates{k2}, insurance_policy{k3},  point_solution(j,7), point_solution(j,8), point_solution(j,9));  
       fprintf(fid,'\n%9s%10s%14s%17s%10.2f%8.2f%8.2f', point_variables{j}, crop_type{k1}, planting_dates{k2}, insurance_policy{k3},  point_solution(j,7), point_solution(j,8), point_solution(j,9));  
    end  
end
fprintf('\n\n');

% display Graph "El Nino Climate Type. Expected Profit vs. CVaR Risk@90%"
clear x;
clear y;
for i=1:1:3
    [x(i)] = functionvalue('cvar_risk', 0.9, matrix_scen_el, [], [], point_solution(1:n_variables_principal,i));
    [y(i)] = functionvalue('linear', [], matrix_objective_el, [], [], point_solution(1:n_variables_principal,i));
end
S = [70; 30; 30];
Col = [[1 0 0]; [0 0 1]; [0 1 0]];

figure
scatter(x,y,S,Col,'filled');
grid on;
xlabel('CVaR Risk with Confidence Level 90% (USD)');
ylabel('Expected Profit (USD)');
text(x(1),y(1),'optimized with CVaR@90%', 'HorizontalAlignment','left', 'VerticalAlignment','top');
text(x(2),y(2),'optimized with PrPen@600', 'HorizontalAlignment','left', 'VerticalAlignment','bottom');
text(x(3),y(3),'optimized with VaR@90%', 'HorizontalAlignment','right', 'VerticalAlignment','top');
title('El Nino Climate Type. Expected Profit vs. CVaR Risk@90%');

% display Graph "El Nino Climate Type. Expected Profit vs. PrPen@600"
clear x;
clear y;
for i=1:1:3
    [x(i)] = functionvalue('pr_pen', 600, matrix_scen_el, [], [], point_solution(1:n_variables_principal,i));
    [y(i)] = functionvalue('linear', [], matrix_objective_el, [], [], point_solution(1:n_variables_principal,i));
end
S = [30; 70; 30];
Col = [[0 0 1]; [1 0 0]; [0 1 0]];

figure
scatter(x,y,S,Col,'filled');
grid on;
xlabel('Probability Exceeding with Threshold = 600 USD (fractions)');
ylabel('Expected Profit (USD)');
text(x(1),y(1),'optimized with CVaR@90%', 'HorizontalAlignment','right', 'VerticalAlignment','top');
text(x(2),y(2),'optimized with PrPen@600', 'HorizontalAlignment','left', 'VerticalAlignment','bottom');
text(x(3),y(3),'optimized with VaR@90%', 'HorizontalAlignment','right', 'VerticalAlignment','top');
title('El Nino Climate Type. Expected Profit vs. PrPen@600');

% display Graph "El Nino Climate Type. Expected Profit vs. VaR Risk@90%"
clear x;
clear y;
for i=1:1:3
    [x(i)] = functionvalue('var_risk', 0.9, matrix_scen_el, [], [], point_solution(1:n_variables_principal,i));
    [y(i)] = functionvalue('linear', [], matrix_objective_el, [], [], point_solution(1:n_variables_principal,i));
end
S = [30; 30; 70];
Col = [[0 1 0]; [0 0 1]; [1 0 0]];

figure
scatter(x,y,S,Col,'filled');
grid on;
xlabel('VaR Risk with Confidence Level 90% (USD)');
ylabel('Expected Profit (USD)');
text(x(1),y(1),'optimized with CVaR@90%', 'HorizontalAlignment','right', 'VerticalAlignment','top');
text(x(2),y(2),'optimized with PrPen@600', 'HorizontalAlignment','left', 'VerticalAlignment','bottom');
text(x(3),y(3),'optimized with VaR@90%', 'HorizontalAlignment','right', 'VerticalAlignment','top');
title('El Nino Climate Type. Expected Profit vs. VaR Risk@90%');


% display Graph "La Nino Climate Type. Expected Profit vs. CVaR Risk@90%"
clear x;
clear y;
for i=1:1:3
    [x(i)] = functionvalue('cvar_risk', 0.9, matrix_scen_la, [], [], point_solution(1:n_variables_principal,i+3));
    [y(i)] = functionvalue('linear', [], matrix_objective_la, [], [], point_solution(1:n_variables_principal,i+3));
end
S = [70; 30; 30];
Col = [[1 0 0]; [0 0 1]; [0 1 0]];

figure
scatter(x,y,S,Col,'filled');
grid on;
xlabel('CVaR Risk with Confidence Level 90% (USD)');
ylabel('Expected Profit (USD)');
text(x(1),y(1),'optimized with CVaR@90%', 'HorizontalAlignment','left', 'VerticalAlignment','top');
text(x(2),y(2),'optimized with PrPen@600', 'HorizontalAlignment','right', 'VerticalAlignment','bottom');
text(x(3),y(3),'optimized with VaR@90%', 'HorizontalAlignment','left', 'VerticalAlignment','top');
title('La Nino Climate Type. Expected Profit vs. CVaR Risk@90%');

% display Graph "La Nino Climate Type. Expected Profit vs. PrPen@600"
clear x;
clear y;
for i=1:1:3
    [x(i)] = functionvalue('pr_pen', 600, matrix_scen_la, [], [], point_solution(1:n_variables_principal,i+3));
    [y(i)] = functionvalue('linear', [], matrix_objective_la, [], [], point_solution(1:n_variables_principal,i+3));
end
S = [30; 70; 30];
Col = [[0 0 1]; [1 0 0]; [0 1 0]];

figure
scatter(x,y,S,Col,'filled');
grid on;
xlabel('Probability Exceeding with Threshold = 600 USD (fractions)');
ylabel('Expected Profit (USD)');
text(x(1),y(1),'optimized with CVaR@90%', 'HorizontalAlignment','left', 'VerticalAlignment','top');
text(x(2),y(2),'optimized with PrPen@600', 'HorizontalAlignment','right', 'VerticalAlignment','bottom');
text(x(3),y(3),'optimized with VaR@90%', 'HorizontalAlignment','left', 'VerticalAlignment','top');
title('La Nino Climate Type. Expected Profit vs. PrPen@600');

% display Graph "La Nino Climate Type. Expected Profit vs. VaR Risk@90%"
clear x;
clear y;
for i=1:1:3
    [x(i)] = functionvalue('var_risk', 0.9, matrix_scen_la, [], [], point_solution(1:n_variables_principal,i+3));
    [y(i)] = functionvalue('linear', [], matrix_objective_la, [], [], point_solution(1:n_variables_principal,i+3));
end
S = [30; 30; 70];
Col = [[0 1 0]; [0 0 1]; [1 0 0]];

figure
scatter(x,y,S,Col,'filled');
grid on;
xlabel('VaR Risk with Confidence Level 90% (USD)');
ylabel('Expected Profit (USD)');
text(x(1),y(1),'optimized with CVaR@90%', 'HorizontalAlignment','left', 'VerticalAlignment','top');
text(x(2),y(2),'optimized with PrPen@600', 'HorizontalAlignment','right', 'VerticalAlignment','bottom');
text(x(3),y(3),'optimized with VaR@90%', 'HorizontalAlignment','left', 'VerticalAlignment','top');
title('La Nino Climate Type. Expected Profit vs. VaR Risk@90%');

% display Graph "Neutral Climate Type. Expected Profit vs. CVaR Risk@90%"
clear x;
clear y;
for i=1:1:3
    [x(i)] = functionvalue('cvar_risk', 0.9, matrix_scen_nt, [], [], point_solution(1:n_variables_principal,i+6));
    [y(i)] = functionvalue('linear', [], matrix_objective_nt, [], [], point_solution(1:n_variables_principal,i+6));
end
S = [70; 30; 30];
Col = [[1 0 0]; [0 0 1]; [0 1 0]];

figure
scatter(x,y,S,Col,'filled');
grid on;
xlabel('CVaR Risk with Confidence Level 90% (USD)');
ylabel('Expected Profit (USD)');
text(x(1),y(1),'optimized with CVaR@90%', 'HorizontalAlignment','left', 'VerticalAlignment','top');
text(x(2),y(2),'optimized with PrPen@600', 'HorizontalAlignment','right', 'VerticalAlignment','bottom');
text(x(3),y(3),'optimized with VaR@90%', 'HorizontalAlignment','left', 'VerticalAlignment','top');
title('Neutral Climate Type. Expected Profit vs. CVaR Risk@90%');

% display Graph "Neutral Climate Type. Expected Profit vs. PrPen@600"
clear x;
clear y;
for i=1:1:3
    [x(i)] = functionvalue('pr_pen', 600, matrix_scen_nt, [], [], point_solution(1:n_variables_principal,i+6));
    [y(i)] = functionvalue('linear', [], matrix_objective_nt, [], [], point_solution(1:n_variables_principal,i+6));
end
S = [30; 70; 30];
Col = [[0 0 1]; [1 0 0]; [0 1 0]];

figure
scatter(x,y,S,Col,'filled');
grid on;
xlabel('Probability Exceeding with Threshold = 600 USD (fractions)');
ylabel('Expected Profit (USD)');
text(x(1),y(1),'optimized with CVaR@90%', 'HorizontalAlignment','left', 'VerticalAlignment','top');
text(x(2),y(2),'optimized with PrPen@600', 'HorizontalAlignment','right', 'VerticalAlignment','bottom');
text(x(3),y(3),'optimized with VaR@90%', 'HorizontalAlignment','left', 'VerticalAlignment','top');
title('Neutral Climate Type. Expected Profit vs. PrPen@600');

% display Graph "Neutral Climate Type. Expected Profit vs. VaR Risk@90%"
clear x;
clear y;
for i=1:1:3
    [x(i)] = functionvalue('var_risk', 0.9, matrix_scen_nt, [], [], point_solution(1:n_variables_principal,i+6));
    [y(i)] = functionvalue('linear', [], matrix_objective_nt, [], [], point_solution(1:n_variables_principal,i+6));
end
S = [30; 30; 70];
Col = [[0 1 0]; [0 0 1]; [1 0 0]];

figure
scatter(x,y,S,Col,'filled');
grid on;
xlabel('VaR Risk with Confidence Level 90% (USD)');
ylabel('Expected Profit (USD)');
text(x(1),y(1),'optimized with CVaR@90%', 'HorizontalAlignment','left', 'VerticalAlignment','top');
text(x(2),y(2),'optimized with PrPen@600', 'HorizontalAlignment','right', 'VerticalAlignment','bottom');
text(x(3),y(3),'optimized with VaR@90%', 'HorizontalAlignment','left', 'VerticalAlignment','top');
title('Neutral Climate Type. Expected Profit vs. VaR Risk@90%');

% display Graph "El Nino Climate Type. Contributions to CVaR Risk@90% of Optimal Solution Components "
clear x;
clear y;
[Increments] = functionincrement('cvar_risk', 0.9, matrix_objective_el, [], [], point_solution(1:n_variables_principal,1));
    j = 0;
    for i=1:1:n_variables_principal
        if abs(Increments(i)) >0
          j = j+1;
          y(j) = Increments(i);
          x(j) = point_variables(i);
        end
    end

  figure
  bar(y,0.1,'g');
  set(gca,'XTick',1:1:length(x));
hAxes = gca;
set(hAxes, 'XTickLabel', x);
grid on;
xlabel('Components of Optimal Solution (Nonzero Conributors)');
ylabel('Contributions to CVaR Risk@90% (USD)');
title('El Nino Climate Type. Contributions to CVaR Risk@90% of Optimal Solution Components');

% display Graph "La Nino Climate Type. Contributions to CVaR Risk@90% of Optimal Solution Components "
clear x;
clear y;
clear Increments;
[Increments] = functionincrement('cvar_risk', 0.9, matrix_objective_la, [], [], point_solution(1:n_variables_principal,4));
    j = 0;
    for i=1:1:n_variables_principal
        if abs(Increments(i)) >0
          j = j+1;
          y(j) = Increments(i);
          x(j) = point_variables(i);
        end
    end

  figure
  bar(y,0.1,'c');
  set(gca,'XTick',1:1:length(x));
hAxes = gca;
set(hAxes, 'XTickLabel', x);
grid on;
xlabel('Components of Optimal Solution (Nonzero Conributors)');
ylabel('Contributions to CVaR Risk@90% (USD)');
title('La Nino Climate Type. Contributions to CVaR Risk@90% of Optimal Solution Components');


% display Graph "Neutral Climate Type. Contributions to CVaR Risk@90% of Optimal Solution Components "
clear x;
clear y;
clear Increments;
[Increments] = functionincrement('cvar_risk', 0.9, matrix_objective_nt, [], [], point_solution(1:n_variables_principal,7));
    j = 0;
    for i=1:1:n_variables_principal
        if abs(Increments(i)) >0
          j = j+1;
          y(j) = Increments(i);
          x(j) = point_variables(i);
        end
    end

  figure
  bar(y,0.1,'r');
  set(gca,'XTick',1:1:length(x));
hAxes = gca;
set(hAxes, 'XTickLabel', x);
grid on;
xlabel('Components of Optimal Solution (Nonzero Conributors)');
ylabel('Contributions to CVaR Risk@90% (USD)');
title('Neutral Climate Type. Contributions to CVaR Risk@90% of Optimal Solution Components');

% display Graph "El Nino Climate Type. Nonzero Components of Optimal Solution Obtained with CVaR Risk@90%
clear x;
clear y;
i = 0;
for j=1:1:n_variables_principal
    if point_solution(j,1)> 0
        i = i+1;
        x(i) = point_variables(j);
        y(i) = point_solution(j,1);
    end  
end
figure
 bar(y,0.1,'g');
  set(gca,'XTick',1:1:length(x));
hAxes = gca;
set(hAxes, 'XTickLabel', x);
grid on;
xlabel('Nonzero Components of Optimal Solution');
ylabel('Planting Area (in acres)');
title('El Nino Climate Type. Nonzero Components of Optimal Solution Obtained with CVaR Risk@90%');

% display Graph "El Nino Climate Type. Nonzero Components of Optimal
% Solution Obtained with PrPen@600
clear x;
clear y;
i = 0;
for j=1:1:n_variables_principal
    if point_solution(j,2)> 0
        i = i+1;
        x(i) = point_variables(j);
        y(i) = point_solution(j,2);
    end  
end
figure
 bar(y,0.1,'b');
  set(gca,'XTick',1:1:length(x));
hAxes = gca;
set(hAxes, 'XTickLabel', x);
grid on;
xlabel('Nonzero Components of Optimal Solution');
ylabel('Planting Area (in acres)');
title('El Nino Climate Type. Nonzero Components of Optimal Solution Obtained with PrPen@600');

% display Graph "El Nino Climate Type. Nonzero Components of Optimal
% Solution Obtained with VaR@90%
clear x;
clear y;
i = 0;
for j=1:1:n_variables_principal
    if point_solution(j,3)> 0
        i = i+1;
        x(i) = point_variables(j);
        y(i) = point_solution(j,3);
    end  
end
figure
 bar(y,0.1,'c');
  set(gca,'XTick',1:1:length(x));
hAxes = gca;
set(hAxes, 'XTickLabel', x);
grid on;
xlabel('Nonzero Components of Optimal Solution');
ylabel('Planting Area (in acres)');
title('El Nino Climate Type. Nonzero Components of Optimal Solution Obtained with VaR@90%');

% display Graph "La Nino Climate Type. Nonzero Components of Optimal Solution Obtained with CVaR Risk@90%

clear x;
clear y;
i = 0;
for j=1:1:n_variables_principal
    if point_solution(j,4)> 0
        i = i+1;
        x(i) = point_variables(j);
        y(i) = point_solution(j,4);
    end  
end
figure
 bar(y,0.1,'m');
  set(gca,'XTick',1:1:length(x));
hAxes = gca;
set(hAxes, 'XTickLabel', x);
grid on;
xlabel('Nonzero Components of Optimal Solution');
ylabel('Planting Area (in acres)');
title('La Nino Climate Type. Nonzero Components of Optimal Solution Obtained with CVaR Risk@90%');


% display Graph "La Nino Climate Type. Nonzero Components of Optimal
% Solution Obtained with PrPen@600

clear x;
clear y;
i = 0;
for j=1:1:n_variables_principal
    if point_solution(j,5)> 0
        i = i+1;
        x(i) = point_variables(j);
        y(i) = point_solution(j,5);
    end  
end
figure
 bar(y,0.1,'y');
  set(gca,'XTick',1:1:length(x));
hAxes = gca;
set(hAxes, 'XTickLabel', x);
grid on;
xlabel('Nonzero Components of Optimal Solution');
ylabel('Planting Area (in acres)');
title('La Nino Climate Type. Nonzero Components of Optimal Solution Obtained with PrPen@600');


% display Graph "La Nino Climate Type. Nonzero Components of Optimal
% Solution Obtained with VaR@90%

clear x;
clear y;
i = 0;
for j=1:1:n_variables_principal
    if point_solution(j,6)> 0
        i = i+1;
        x(i) = point_variables(j);
        y(i) = point_solution(j,6);
    end  
end
figure
 bar(y,0.1,'r');
  set(gca,'XTick',1:1:length(x));
hAxes = gca;
set(hAxes, 'XTickLabel', x);
grid on;
xlabel('Nonzero Components of Optimal Solution');
ylabel('Planting Area (in acres)');
title('La Nino Climate Type. Nonzero Components of Optimal Solution Obtained with VaR@90%');


% display Graph "Neutral Climate Type. Nonzero Components of Optimal Solution Obtained with CVaR Risk@90%

clear x;
clear y;
i = 0;
for j=1:1:n_variables_principal
    if point_solution(j,7)> 0
        i = i+1;
        x(i) = point_variables(j);
        y(i) = point_solution(j,7);
    end  
end
figure
 bar(y,0.1,'g');
  set(gca,'XTick',1:1:length(x));
hAxes = gca;
set(hAxes, 'XTickLabel', x);
grid on;
xlabel('Nonzero Components of Optimal Solution');
ylabel('Planting Area (in acres)');
title('Neutral Climate Type. Nonzero Components of Optimal Solution Obtained with CVaR Risk@90%');


% display Graph "Neutral Climate Type. Nonzero Components of Optimal
% Solution Obtained with PrPen@600

clear x;
clear y;
i = 0;
for j=1:1:n_variables_principal
    if point_solution(j,8)> 0
        i = i+1;
        x(i) = point_variables(j);
        y(i) = point_solution(j,8);
    end  
end
figure
 bar(y,0.1,'b');
  set(gca,'XTick',1:1:length(x));
hAxes = gca;
set(hAxes, 'XTickLabel', x);
grid on;
xlabel('Nonzero Components of Optimal Solution');
ylabel('Planting Area (in acres)');
title('Neutral Climate Type. Nonzero Components of Optimal Solution Obtained with PrPen@600');


% display Graph "Neutral Climate Type. Nonzero Components of Optimal
% Solution Obtained with VaR@90%

clear x;
clear y;
i = 0;
for j=1:1:n_variables_principal
    if point_solution(j,9)> 0
        i = i+1;
        x(i) = point_variables(j);
        y(i) = point_solution(j,9);
    end  
end
figure
 bar(y,0.1,'c');
  set(gca,'XTick',1:1:length(x));
hAxes = gca;
set(hAxes, 'XTickLabel', x);
grid on;
xlabel('Nonzero Components of Optimal Solution');
ylabel('Planting Area (in acres)');
title('Neutral Climate Type. Nonzero Components of Optimal Solution Obtained with VaR@90%');

%=======================================================================|
%American Optimal Decisions, Inc. Copyright                             |
%Copyright ©American Optimal Decisions, Inc. 2007-2015.                 |
%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.              |
%=======================================================================|