It turns out that calling the network solver to compute a shortest path in a directed acyclic network is a bit of overkill. You can instead apply dynamic programming as follows, and this approach uses less memory:
proc optmodel;
/* read data */
set OBS;
str obs_id {OBS};
num pd {OBS};
num def {OBS};
num cust_cnt {OBS};
read data have into OBS=[_N_] obs_id=customer pd def cust_cnt;
/* build network */
num dummy_obs init card(OBS)+1;
OBS = OBS union {dummy_obs};
num num_cust_total = sum {i in OBS diff {dummy_obs}} cust_cnt[i];
num num_cust_lb = ceil( &pct_lb*num_cust_total);
num num_cust_ub = floor(&pct_ub*num_cust_total);
set <num,num> OBS_PAIRS init {};
num num_cust_ij {OBS_PAIRS};
num count;
for {i in OBS} do;
count = 0;
for {j in OBS: i < j} do;
count = count + cust_cnt[j-1];
if count < num_cust_lb then continue;
else if count > num_cust_ub then leave;
else do;
OBS_PAIRS = OBS_PAIRS union {<i,j>};
num_cust_ij[i,j] = count;
end;
end;
end;
/* avg_def[i,j] = (sum {k in i..j-1} def[k] * cust_cnt[k]) / num_cust_ij[i,j] */
num avg_def {<i,j> in OBS_PAIRS} = (
if <i,j-1> not in OBS_PAIRS then (sum {k in i..j-1} def[k] * cust_cnt[k]) / num_cust_ij[i,j]
else (avg_def[i,j-1] * num_cust_ij[i,j-1] + def[j-1] * cust_cnt[j-1]) / num_cust_ij[i,j]
);
num error {<i,j> in OBS_PAIRS} = sum {k in i..j-1} abs(pd[k] - avg_def[i,j]);
set OBS_BUCKETS = OBS cross 1..&num_buckets+1;
/* forward pass to compute shortest-path distances */
num dist_from_source {OBS_BUCKETS} init constant('BIG');
dist_from_source[1,1] = 0;
for {b in 1..&num_buckets} do;
put b=;
if b = 1 then do;
for {<(1),j> in OBS_PAIRS} do;
if dist_from_source[j,b+1] > dist_from_source[1,b] + error[1,j] then do;
dist_from_source[j,b+1] = dist_from_source[1,b] + error[1,j];
end;
end;
end;
else if b = &num_buckets then do;
for {<i,(dummy_obs)> in OBS_PAIRS} do;
if dist_from_source[dummy_obs,b+1] > dist_from_source[i,b] + error[i,dummy_obs] then do;
dist_from_source[dummy_obs,b+1] = dist_from_source[i,b] + error[i,dummy_obs];
end;
end;
end;
else do;
for {<i,j> in OBS_PAIRS: i ne 1 and j ne dummy_obs} do;
if dist_from_source[j,b+1] > dist_from_source[i,b] + error[i,j] then do;
dist_from_source[j,b+1] = dist_from_source[i,b] + error[i,j];
end;
end;
end;
end;
/* backward pass to extract shortest path */
set <num,num,num> PATH init {}; /* order, from, to */
num curr_i init dummy_obs;
num curr_j;
for {b in &num_buckets...1 by -1} do;
curr_j = curr_i;
for {<i,(curr_j)> in OBS_PAIRS} do;
if dist_from_source[curr_j,b+1] = dist_from_source[i,b] + error[i,curr_j] then do;
curr_i = i;
leave;
end;
end;
PATH = {<b,curr_i,curr_j>} union PATH;
end;
put dist_from_source[dummy_obs,&num_buckets+1]=;
put PATH;
/* highest pd value in each bucket */
num cutoff {1..&num_buckets};
num avg_def_b {1..&num_buckets};
num error_b {1..&num_buckets};
num num_cust_b {1..&num_buckets};
num pct_cust_b {b in 1..&num_buckets} = num_cust_b[b] / num_cust_total;
set <str> OBS_b {1..&num_buckets};
for {<order,from,to> in PATH} do;
cutoff[order] = pd[to-1];
avg_def_b[order] = avg_def[from,to];
error_b[order] = error[from,to];
num_cust_b[order] = num_cust_ij[from,to];
OBS_b[order] = setof {k in from..to-1} obs_id[k];
end;
print dist_from_source[dummy_obs,&num_buckets+1] cutoff avg_def_b error_b num_cust_b pct_cust_b;
/* observations in each bucket */
/* for {b in 1..&num_buckets} put OBS_b[b]=;*/
quit;
