Does the following PROC OPTMODEL code do what you want?
proc optmodel;
/* declare parameters and read data */
set GROUPS = /Test Control/;
num groupSize {GROUPS} = [18, 5];
set OBS;
str PatientGroup {OBS};
num Balance {OBS};
num Risk {OBS};
str prod {OBS};
read data indata into OBS=[_N_] PatientGroup Balance Risk prod;
set PATIENT_GROUPS = setof {i in OBS} PatientGroup[i];
set PRODS = setof {i in OBS} prod[i];
set OBS_p {PATIENT_GROUPS} init {};
num patientCount {PATIENT_GROUPS} init 0;
num sumBalance {PATIENT_GROUPS} init 0;
num sumRisk {PATIENT_GROUPS} init 0;
num sumProd {PATIENT_GROUPS, PRODS} init 0;
str p_this, pr_this;
for {i in OBS} do;
p_this = PatientGroup[i];
OBS_p[p_this] = OBS_p[p_this] union {i};
patientCount[p_this] = patientCount[p_this] + 1;
sumBalance[p_this] = sumBalance[p_this] + Balance[i];
sumRisk[p_this] = sumRisk[p_this] + Risk[i];
pr_this = prod[i];
sumProd[p_this, pr_this] = sumProd[p_this, pr_this] + 1;
end;
num totalBalance = sum {p in PATIENT_GROUPS} sumBalance[p];
num totalRisk = sum {p in PATIENT_GROUPS} sumRisk[p];
num totalProd {pr in PRODS} = sum {p in PATIENT_GROUPS} sumProd[p,pr];
/* Assign[p,g] = 1 if patient group p is assigned to group g, 0 otherwise */
var Assign {PATIENT_GROUPS, GROUPS} binary;
/* assign each patient group to exactly one group */
con AssignOnce {p in PATIENT_GROUPS}:
sum {g in GROUPS} Assign[p,g] = 1;
/* assign groupSize[g] patient groups to group g */
con GroupSizeCon {g in GROUPS}:
sum {p in PATIENT_GROUPS} Assign[p,g] = groupSize[g];
/* penalize absolute difference from fair share for Balance */
var BalanceErrorPlus {GROUPS} >= 0;
var BalanceErrorMinus {GROUPS} >= 0;
con BalanceErrorCon {g in GROUPS}:
(sum {p in PATIENT_GROUPS} sumBalance[p] * Assign[p,g]) / totalBalance - groupSize[g] / card(PATIENT_GROUPS)
= BalanceErrorPlus[g] - BalanceErrorMinus[g];
impvar BalanceError = sum {g in GROUPS} (BalanceErrorPlus[g] + BalanceErrorMinus[g]);
/* penalize absolute difference from fair share for Risk */
var RiskErrorPlus {GROUPS} >= 0;
var RiskErrorMinus {GROUPS} >= 0;
con RiskErrorCon {g in GROUPS}:
(sum {p in PATIENT_GROUPS} sumRisk[p] * Assign[p,g]) / totalRisk - groupSize[g] / card(PATIENT_GROUPS)
= RiskErrorPlus[g] - RiskErrorMinus[g];
impvar RiskError = sum {g in GROUPS} (RiskErrorPlus[g] + RiskErrorMinus[g]);
/* penalize absolute difference from fair share for Prod */
var ProdErrorPlus {GROUPS, PRODS} >= 0;
var ProdErrorMinus {GROUPS, PRODS} >= 0;
con ProdErrorCon {g in GROUPS, pr in PRODS}:
(sum {p in PATIENT_GROUPS} sumProd[p,pr] * Assign[p,g]) / totalProd[pr] - groupSize[g] / card(PATIENT_GROUPS)
= ProdErrorPlus[g,pr] - ProdErrorMinus[g,pr];
impvar ProdError = sum {g in GROUPS, pr in PRODS} (ProdErrorPlus[g,pr] + ProdErrorMinus[g,pr]);
/* minimize total error */
min TotalError = BalanceError + RiskError + ProdError;
/* call MILP solver */
solve;
/* create output data */
impvar MeanBalance {g in GROUPS} =
(sum {p in PATIENT_GROUPS} sumBalance[p] * Assign[p,g]) / (sum {p in PATIENT_GROUPS} patientCount[p] * Assign[p,g]);
impvar MeanRisk {g in GROUPS} =
(sum {p in PATIENT_GROUPS} sumRisk[p] * Assign[p,g]) / (sum {p in PATIENT_GROUPS} patientCount[p] * Assign[p,g]);
impvar MeanProd {g in GROUPS, pr in PRODS} =
(sum {p in PATIENT_GROUPS} sumProd[p,pr] * Assign[p,g]) / (sum {p in PATIENT_GROUPS} patientCount[p] * Assign[p,g]);
create data GroupOutdata from [Group] groupSize MeanBalance MeanRisk
{pr in PRODS} <col(pr)=MeanProd[Group,pr]>;
str assignedGroup {PATIENT_GROUPS};
for {p in PATIENT_GROUPS} do;
for {g in GROUPS: Assign[p,g].sol > 0.5} do;
assignedGroup[p] = g;
leave;
end;
end;
create data PatientGroupOutdata from [PatientGroup=p] assignedGroup patientCount sumBalance sumRisk
{pr in PRODS} <col(pr)=sumProd[p,pr]>;
quit;
