Open cast mining problem - Linear Matrix Inequality method

Problem statement link: 

 http://support.sas.com/documentation/cdl/en/ormpex/67518/HTML/default/viewer.htm#ormpex_ex14_sect001.htm

Mining problem is a famous one in operations research, in this post I show my solution for this problem using linear matrix inequality representation as follows:

where f = value - cost for block i and n = 30. Consider the constraints to extract block #16,

then the sum of these constraints will be as follows, 

rearranging this term,

this constraint is the first row of the product AX i.e,,

using the same way to describe the constraints for the other blocks, then this ends up with the above formulation. 

Results:  Optimal value = 17500,  Blocks to be extracted: 0, 1 ,2 , 4, 5, 6, 8, 9, 10, 16, 17, 19, 20, 25.

Heuristic solution by Matlab:

%solving without solvers - heuristic solution

clear
clc


cost = [3000*ones(1,16), 6000*ones(1,9), 8000*ones(1,4), 10000];

value = 2000* [1.5, 1.5, 1.5, 0.75,...
    1.5, 2, 1.5, 0.75,...
    1, 1, 0.75, 0.5,...
    0.75, 0.75, 0.5, 0.25,...
    4, 4, 2,...
    3, 3, 1,...
    2, 2, 0.5,...
    12, 6,...
    5, 4, ...
    6];

neighbors_16=[[16,0], [16,1], [16,4], [16,5]];
neighbors_17=[[17,1], [17,2], [17,5], [17,6]];
neighbors_25=[[25,16], [25,17], [25,19], [25,20],...
    [25,0], [25,1], [25,2], [25,4], [25,5], [25,6],...
    [25,8], [25,9], [25,10]];

neighbors_26=[[26,17], [26,18], [26,20], [26,21],...
    [26,1], [26,2], [26,3], [26,5], [26,6], [26,7],...
    [26,9], [26,10], [26,11]];

neighbors_27=[[27,19], [27,20], [27,22], [27,23],...
    [27,4], [27,5], [27,6], [27,8], [27,9], [27,10],...
    [27,12], [27,13], [27,14]];

neighbors_28=[[28,20], [28,21], [28,23], [27,24],...
      [27,5], [28,6], [28,7], [28,9], [28,10], [28,11],...
      [28,13], [27,14], [27,15]];

neighbors_29=[[29, 0], [29, 1], [29,2], [29,3],...
    [29, 4], [29, 5], [29,6], [29,7],...
    [29, 8], [29, 9], [29,10], [29,11],...
    [29, 12], [29, 13], [29,14], [29,15],...
    [29, 16], [29, 17], [29,18], [29,19],...
    [29, 20], [29, 21], [29,22], [29,23], [29,24]];


edges_total=[neighbors_16 zeros(1,50-length(neighbors_16));...
   neighbors_17 zeros(1,50-length(neighbors_17));...
   neighbors_25 zeros(1,50-length(neighbors_25));...
   neighbors_26 zeros(1,50-length(neighbors_26));...
   neighbors_27 zeros(1,50-length(neighbors_27));...
   neighbors_28 zeros(1,50-length(neighbors_28));...
   neighbors_29 zeros(1,50-length(neighbors_29))];
    

p_total=value-cost;

[p_max_total,l_max_total]=max(p_total);

index=find(value-cost > 0);

block_index=index-1;

value=value(index);

cost=cost(index);

p=value-cost; %non zero profits

[p_max,l_max]=max(p);

edges = edges_total(l_max-1,[1:length(neighbors_25)]);

x=edges(2:2:(length(edges)));

optimal_obj_value=sum(p_total(x+1))+max(p)

blocks_extracted=sort([x (l_max_total-1)])