About RobPratt

RobPratt · ‎04-16-2017

You can use the EXPAND statement to see how the objective function you defined gets populated: Var x >= 0 Var y >= 0 Minimize MinWgtDis=SQRT((-x + 100)**2 + (-y + 100)**2)*425 + SQRT((-x + 80)**2 + (-y + 80) **2)*425 + SQRT((-x + 86)**2 + (-y + 86)**2)*320 + SQRT((-x + 40)**2 + (-y + 40)**2)*320 + SQRT((-x + 20)**2 + (-y + 20)**2)*220 + SQRT((-x + 60)**2 + (-y + 60)**2)*220 You can see that the formula for MinWgtDis is not really what you intended. Here's a modification of your code that instead does what you want: Proc optmodel; set city = {'BO','PR','SP'}; set Geo = {'Xi','Yi'}; number Cor{city,Geo} =[100 80 86 40 20 60]; number Wgt{city}=[425 320 220]; var x {Geo} >= 0 ; Minimize MinWgtDis = sum{i in city} sqrt(sum {j in Geo} (cor[i,j] - x[j])^2)* Wgt[i]; solve; print MinWgtDis x; quit;

RobPratt · ‎04-16-2017

Here's one way to solve it as an unconstrained problem with nonlinear objective and two variables: data indata; input city $12. Wgt X Y; datalines; Boston 425 100 80 Providence 320 86 40 Springfield 220 20 60 ; proc sgplot data=indata; bubble x=x y=y size=wgt / datalabel=city; run; proc optmodel; set <str> CITIES; num wgt {CITIES}; num xc {CITIES}; num yc {CITIES}; read data indata into CITIES=[city] wgt xc=x yc=y; var X, Y; min MinWgtDis = sum {i in CITIES} (sqrt((xc[i] - X)^2 + (yc[i] - Y)^2)) * wgt[i]; solve; print MinWgtDis X Y; create data sganno from function='text' drawspace='datavalue' label='facility' x1=X y1=Y; quit; proc sgplot data=indata sganno=sganno; bubble x=x y=y size=wgt / datalabel=city; run; The resulting plot provides a sanity check on the optimal solution returned.

RobPratt · ‎04-16-2017

Because PROC PHREG is in SAS/STAT, please post in the stat procs community instead.

RobPratt · ‎04-11-2017

There is a hard limit of 2^31 - 1 = 2,147,483,647 variables, constraints, or nonzero coefficients. Other than that, the only limit depends on the amount of memory you have provided to SAS.

RobPratt · ‎04-10-2017

You can formulate and solve the problem with PROC OPTMODEL in SAS/OR: proc optmodel; set ISET = 1..3; set JSET = 1..30000; num p {ISET, JSET} = rand('UNIFORM'); num b {ISET} = card(JSET)*rand('UNIFORM'); var X {ISET, JSET} binary; max Z = sum {i in ISET, j in JSET} p[i,j] * X[i,j]; con Con1 {j in JSET}: sum {i in ISET} X[i,j] = 1; con Con2 {i in ISET}: sum {j in JSET} X[i,j] <= b[i]; /* call (default) MILP solver */ solve; /* call LP solver with (default) dual simplex algorithm */ solve with lp relaxint; /* call LP solver with network simplex algorithm */ solve with lp relaxint / algorithm=ns; quit; Because the problem has a pure network structure, the third solve performs the best, solving in less than a second. Due to total unimodularity of the resulting constraint matrix, you can relax the integrality of x and an optimal solution will automatically take integer values. You can also model the problem in terms of a bipartite directed network, where the left side consists of the i nodes, each with a supply of at most b[i], and the right side consists of the j nodes, each with a demand of exactly 1. The (binary) variable x[i,j] represents the flow from node i to node j. You can negate the link weights (because you want to maximize) and call the minimum-cost network flow algorithm either via the network solver in PROC OPTMODEL or via PROC OPTNET, also in SAS/OR.

RobPratt · ‎04-10-2017

Here's how you can solve it with the NLP solver in PROC OPTMODEL: proc optmodel; var x >= -10 <= 10 init 1; var y >= -6 <= 6 init 1; min F = 0.1*x*y; con Mycon: x**2 + y**2 <= (5+2.2*cos(10*atan(x/y)))**2; solve with nlp / ms; print x y; quit;

RobPratt · ‎04-05-2017

The error message indicates that the upper bound of the constraint is missing. You can use the EXPAND statement to see the constraint that results from your declaration: expand Max_Budget; Also, it might be simpler to read the data into PROC OPTMODEL directly from subset_Table instead of calling PROC TRANSPOSE and PROC SQL.

RobPratt · ‎03-17-2017

Glad to help. By the way, the new secondary objective also naturally minimizes the number of clusters. So here is an even simpler and shorter version that calls the solver only once. An optimal solution will have the clusters in nonincreasing order of node weight, with all the empty clusters at the end, so I removed the relabeling section. For your real data, it is still conceivable that the two-objective version solves faster, so please try both. proc optmodel; /* declare parameters and read data */ num max_distance = 70; set NODES; str state {NODES}; str city {NODES}; num latitude {NODES}; num longitude {NODES}; num node_weight {NODES}; read data indata into NODES=[_N_] state city latitude longitude node_weight; num distance {i in NODES, j in NODES: i < j} = geodist(latitude[i], longitude[i], latitude[j], longitude[j], 'DM'); set CLUSTERS init NODES; /* UseCluster[k] = 1 if cluster k is used, 0 otherwise */ var UseCluster {CLUSTERS} binary; /* AssignNodeToCluster[i,k] = 1 if node i is assigned to cluster k, 0 otherwise */ var AssignNodeToCluster {NODES, CLUSTERS} binary; /* assign each node to exactly one cluster */ con AssignOnce {i in NODES}: sum {k in CLUSTERS} AssignNodeToCluster[i,k] = 1; /* if AssignNodeToCluster[i,k] = 1 then UseCluster[k] = 1 */ con AssignToOpenCluster {i in NODES, k in CLUSTERS}: AssignNodeToCluster[i,k] <= UseCluster[k]; /* cannot assign i and j to same cluster k if distance[i,j] > max_distance */ set CONFLICTS = {i in NODES, j in NODES: i < j and distance[i,j] > max_distance}; con Conflict {<i,j> in CONFLICTS, k in CLUSTERS}: AssignNodeToCluster[i,k] + AssignNodeToCluster[j,k] <= UseCluster[k]; /* declare implicit variable to capture total node weight for each cluster */ impvar NodeWeightPerCluster {k in CLUSTERS} = sum {i in NODES} node_weight[i] * AssignNodeToCluster[i,k]; /* objective: move nodes to higher weight clusters if possible */ min MaxNodeWeights = sum {k in CLUSTERS} k * NodeWeightPerCluster[k]; /* call MILP solver */ solve with MILP / decomp=(method=set); num assigned_cluster {NODES}; for {i in NODES} do; for {k in CLUSTERS: AssignNodeToCluster[i,k].sol > 0.5} do; assigned_cluster[i] = k; leave; end; end; /* save solution to SAS data set */ create data outdata from [node] state city latitude longitude node_weight cluster=assigned_cluster; quit; /* print solution by cluster */ proc sort data=outdata; by cluster; run; proc print data=outdata; by cluster; sum node_weight; run; proc sgplot data=outdata; scatter y=latitude x=longitude / datalabel=city group=cluster; run;

RobPratt · ‎03-17-2017

After further thought, here's a better (and simpler) choice for a secondary objective that will now move nodes to higher weight clusters if possible: /* secondary objective: maximize weighted sum of cluster node weights */ max MaxNodeWeights = sum {k in CLUSTERS} k * NodeWeightPerCluster[k]; /* call MILP solver to optimize secondary objective */ solve with MILP / primalin; for {i in NODES} do; for {k in CLUSTERS: AssignNodeToCluster[i,k].sol > 0.5} do; assigned_cluster[i] = k; leave; end; end; In particular, it does assign Charleston to the "right" cluster. Your suggested approach with cliques is interesting. Note that the network algorithms from PROC OPTNET are also available in PROC OPTMODEL via the SOLVE WITH NETWORK statement.

RobPratt · ‎03-16-2017

The code below calls the MILP solver twice, once for each objective. The first solve uses the decomposition algorithm, which exploits the fact that cluster labels are arbitrary and hence the subproblem has identical blocks. (For an amusing application to wedding planning, see this blog post.) The second solve encourages clusters with high total node weights by maximizing the minimum total node weight per cluster. proc optmodel; /* declare parameters and read data */ num max_distance = 70; set NODES; str state {NODES}; str city {NODES}; num latitude {NODES}; num longitude {NODES}; num node_weight {NODES}; read data indata into NODES=[_N_] state city latitude longitude node_weight; num distance {i in NODES, j in NODES: i < j} = geodist(latitude[i], longitude[i], latitude[j], longitude[j], 'DM'); set CLUSTERS init NODES; /* UseCluster[k] = 1 if cluster k is used, 0 otherwise */ var UseCluster {CLUSTERS} binary; /* AssignNodeToCluster[i,k] = 1 if node i is assigned to cluster k, 0 otherwise */ var AssignNodeToCluster {NODES, CLUSTERS} binary; /* primary objective: minimize number of used clusters */ min NumUsedClusters = sum {k in CLUSTERS} UseCluster[k]; /* assign each node to exactly one cluster */ con AssignOnce {i in NODES}: sum {k in CLUSTERS} AssignNodeToCluster[i,k] = 1; /* if AssignNodeToCluster[i,k] = 1 then UseCluster[k] = 1 */ con AssignToOpenCluster {i in NODES, k in CLUSTERS}: AssignNodeToCluster[i,k] <= UseCluster[k]; /* cannot assign i and j to same cluster k if distance[i,j] > max_distance */ set CONFLICTS = {i in NODES, j in NODES: i < j and distance[i,j] > max_distance}; con Conflict {<i,j> in CONFLICTS, k in CLUSTERS}: AssignNodeToCluster[i,k] + AssignNodeToCluster[j,k] <= UseCluster[k]; /* declare implicit variable to capture total node weight for each cluster */ impvar NodeWeightPerCluster {k in CLUSTERS} = sum {i in NODES} node_weight[i] * AssignNodeToCluster[i,k]; /* call MILP solver to optimize primary objective */ solve with MILP / decomp=(method=set); /* relabel clusters from 1 to the minimum number necessary */ num assigned_cluster {NODES}; for {i in NODES} do; for {k in CLUSTERS: AssignNodeToCluster[i,k].sol > 0.5} do; assigned_cluster[i] = k; leave; end; end; num new_cluster {1..card(NODES)}; num new_k init 0; for {k in CLUSTERS: UseCluster[k].sol > 0.5} do; new_k = new_k + 1; new_cluster[k] = new_k; end; for {i in NODES, k in CLUSTERS} AssignNodeToCluster[i,k] = 0; for {i in NODES} do; assigned_cluster[i] = new_cluster[assigned_cluster[i]]; AssignNodeToCluster[i,assigned_cluster[i]] = 1; end; CLUSTERS = 1..new_k; for {k in CLUSTERS} fix UseCluster[k] = 1; /* secondary objective: maximize minimum node weight across clusters */ var MinNodeWeight >= 0 init min {k in CLUSTERS} NodeWeightPerCluster[k].sol; max MaxMin = MinNodeWeight; con MaxMinDef {k in CLUSTERS}: MinNodeWeight <= NodeWeightPerCluster[k]; /* call MILP solver to optimize secondary objective */ solve with MILP / primalin; for {i in NODES} do; for {k in CLUSTERS: AssignNodeToCluster[i,k].sol > 0.5} do; assigned_cluster[i] = k; leave; end; end; /* save solution to SAS data set */ create data outdata from [node] state city latitude longitude node_weight cluster=assigned_cluster; quit; /* print and plot solution by cluster */ proc sort data=outdata; by cluster; run; proc print data=outdata; by cluster; sum node_weight; run; proc sgplot data=outdata; scatter y=latitude x=longitude / datalabel=city group=cluster; run; On my machine running SAS/OR 14.2, I get the following 15 clusters in about 2 seconds: cluster=1 Obs node state city latitude longitude node_weight 1 1 WV Gassaway 38.6732 -80.7748 30 2 2 WV Quinwood 38.0576 -80.7068 40 cluster 70 cluster=2 Obs node state city latitude longitude node_weight 3 3 WV Charleston 38.3498 -81.6326 55 4 4 WV Kenova 38.3990 -82.5782 25 5 5 WV Crum 37.9057 -82.4460 35 cluster 115 cluster=3 Obs node state city latitude longitude node_weight 6 6 AR Fort Smith 35.3859 -94.3985 114 cluster=4 Obs node state city latitude longitude node_weight 7 7 IN Indianapolis 39.7684 -86.1581 28 cluster=5 Obs node state city latitude longitude node_weight 8 8 IN Jasper 38.3914 -86.9311 90 9 9 IN New Albany 38.2856 -85.8241 34 10 10 IN Saint Meinrad 38.1711 -86.8092 67 11 11 IN Scottsburg 38.6856 -85.7702 260 12 12 KY Louisville 38.2527 -85.7585 299 cluster 750 cluster=6 Obs node state city latitude longitude node_weight 13 13 KY Winchester 37.9901 -84.1797 51 cluster=7 Obs node state city latitude longitude node_weight 14 14 NC Charlotte 35.2271 -80.8431 294 cluster=8 Obs node state city latitude longitude node_weight 15 16 NC Horse Shoe 35.343 -82.5567 49 cluster=9 Obs node state city latitude longitude node_weight 16 15 NC High Point 35.9557 -80.0053 76 17 17 NC Kernersville 36.1199 -80.0737 150 18 18 NC Liberty 35.8535 -79.5717 89 cluster 315 cluster=10 Obs node state city latitude longitude node_weight 19 20 NJ Secaucus 40.7895 -74.0565 202 20 21 NJ South River 40.4465 -74.3860 56 cluster 258 cluster=11 Obs node state city latitude longitude node_weight 21 22 NY Geneva 42.868 -76.9856 51 cluster=12 Obs node state city latitude longitude node_weight 22 23 OH Grove City 39.8815 -83.0930 46 23 24 OH Groveport 39.8532 -82.8883 63 cluster 109 cluster=13 Obs node state city latitude longitude node_weight 24 25 OH Youngstown 41.0998 -80.6495 55 cluster=14 Obs node state city latitude longitude node_weight 25 26 PA Corry 41.9203 -79.6403 35 cluster=15 Obs node state city latitude longitude node_weight 26 19 NC Statesville 35.7826 -80.8873 56 27 27 VA Galax 36.6612 -80.9240 43 cluster 99 2393 You can see that Charleston ends up in the "wrong" cluster. The maximin approach essentially ignores all clusters whose total node weights are higher than the minimum, which is 28 in this case. By the way, how many locations are in your real data?

RobPratt · ‎03-15-2017

Please attach it to the discussion.

RobPratt · ‎03-15-2017

1. For help with the FASTCLUS procedure in SAS/STAT, please post in the SAS Statistical Procedures community. 2. The main pitfall of this workaround is that it overconstrains the problem and hence artifically increases the total distance from locations to clusters. An alternative that can capture exactly the constraint you really want is to use the MILP solver in SAS/OR. If you provide data, I'll take a look.

RobPratt · ‎03-07-2017

Those notes arise because the same dummy parameter appears as both a constraint index and a summation index within the constraint. I don't know your optimization model, but I suspect you want the following instead: con RR_BINARY_CONSTRAINT{<ITEMID,STOREID> in Budget_dim}: sum{<(ITEMID),(STOREID),RR> in SKU_STORE_RR_DIM} SKU_STORE_RR_BIN[ITEMID,STOREID,RR]=1; con BUDGET_CONSTRAINT{<ITEMID,STOREID> in Budget_dim}: sum{<(ITEMID),(STOREID),RR> in SKU_STORE_RR_DIM} SKU_STORE_RR_BIN[ITEMID,STOREID,RR]*InvCost[ITEMID,STOREID,RR]*Avg_Onhand[ITEMID,STOREID,RR] <= Budget[ITEMID,STOREID]; At least that will avoid the notes. To see whether the new declarations correctly express your desired model, you can use the EXPAND statement. For more explanation of this implicit slice, see Manpower Planning and other examples in the Mathematical Programming Examples documentation.

RobPratt · ‎02-23-2017

If you have some priorities among the objectives, you can minimize them one at a time sequentially, adding a new optimality constraint at each step. See the SECONDARY OBJECTIVE section of this SAS Global Forum 2014 paper. If you want to perform multiobjective optimization and output a Pareto frontier, you can use PROC OPTLSO. To avoid the failure and warning, you can either set the variables to an initial feasible solution before calling the solver or use the MULTISTART option in the SOLVE WITH NLP statement.

RobPratt · ‎02-23-2017

Your MIN statement declares several objectives (one for each b), but you have not specified which one you want to use. Either declare only one objective or use the OBJ clause in the SOLVE statement to pick one.

Online Status	Offline
Date Last Visited	a week ago