Excessively long time to create models for large systems

[pranav]

Hi,

I am trying to build a linear program to solve economic dispatch problem with huge data set (9000 bus system). Building the model takes a very long time but there are no errors displayed. Even after 2 hours, the model is not built. However the same problem for a smaller system (73 bus) builds the model and solves within 10 minutes. Please let me know if there is any way to speed up the process. I am using Python 3.4.3 and Pyomo 4.0.9638. Thank you.

Regards,

Pranav  

[Gabriel Hackebeil]

The easiest thing to try is running Pyomo through PyPy.

Next, how are you loading your data? Is it through DAT files? If so, there might be faster ways, but you should try profiling before doing anything else. The pyomo command has a few options that should help.

The “--report-timing” option will print the individual construction time for your model components as each is built. If you add this option but don’t see any output for a long time, it likely means this time is being spent in the DAT file parser. It is a good option to see if there are one or two constraints that you need to focus on speeding up.

The “--profile N” option (where N is some number like 100) will show even more profile information. This will likely slow down your code so I would only try it on smaller instances to see if there are any glaring hot spots in your code.

Gabe

--
You received this message because you are subscribed to the Google Groups "Pyomo Forum" group.
To unsubscribe from this group and stop receiving emails from it, send an email to pyomo-forum...@googlegroups.com.
For more options, visit https://groups.google.com/d/optout.

[Cristián Serpell]

If your problem is during formulation, instead of preprocessing or LP writing, we have seen huge improvements in time just modifying the way we write our constraints. Try to avoid innecessary loops, try to get all your "if" comparisons out of the loops, try to use pythonic way to iterate, instead of raw "for"s.

I hope it helps,

Cristián

[pranav]

Hi,

I tried using the "report_timing=True" in the create() method to track the time it takes to build each component of the model.

The result is given below.

70 Declarations: bus branch gen inst ld wind BUSxTIME WINDxTIME GENxTIME BRANCHx                                                   TIME BRANCHxBUS bus_num ld_bus wind_bus branch_fbus branch_tbus linex branch_b b                                                   ranch_maxcap rateC gen_bus gen_pmax gen_pmin gen_heatrate gen_fuelcost gen_sucos                                                   t gen_sdcost gen_nlcost gen_mindowntime gen_minuptime gen_hourlyramp gen_spramp                                                    gen_nspramp gen_susdramp gen_zone1 gen_zone2 gen_zone3 gen_zone4 res_req res_req                                                   _zone1 res_req_zone2 res_req_zone3 res_req_zone4 hourly_load hourly_wind gen_sta                                                   tus gen_vstatus gen_wstatus gen_contingency branch_contingency net_interchange d                                                   emand_multiplier vang ptdf masterU masterV masterW gen_supply sp_res nsp_res inj                                                   ection negativeviolation positiveviolation c_index c obj conlinecap3 conlinecap4                                                    consysbalptdf_rule connodebalance

      0.02 seconds required to construct component=bus; 7726 indicies total

      0.03 seconds required to construct component=branch; 9000 indicies total

         0 seconds required to construct component=gen; 933 indicies total

         0 seconds required to construct component=inst; 24 indicies total

      0.01 seconds required to construct component=ld; 4294 indicies total

         0 seconds required to construct component=wind; 77 indicies total

      0.82 seconds required to construct component=BUSxTIME; 185424 indicies total

      0.01 seconds required to construct component=WINDxTIME; 1848 indicies total

      0.09 seconds required to construct component=GENxTIME; 22392 indicies total

      0.93 seconds required to construct component=BRANCHxTIME; 216000 indicies total

    311.17 seconds required to construct component=BRANCHxBUS; 69534000 indicies total

      0.02 seconds required to construct component=bus_num; 7726 indicies total

      0.01 seconds required to construct component=ld_bus; 4294 indicies total

         0 seconds required to construct component=wind_bus; 77 indicies total

      0.02 seconds required to construct component=branch_fbus; 9000 indicies total

      0.02 seconds required to construct component=branch_tbus; 9000 indicies total

      0.02 seconds required to construct component=linex; 9000 indicies total

      0.02 seconds required to construct component=branch_b; 9000 indicies total

      0.02 seconds required to construct component=branch_maxcap; 9000 indicies total

      0.02 seconds required to construct component=rateC; 9000 indicies total

         0 seconds required to construct component=gen_bus; 933 indicies total

         0 seconds required to construct component=gen_pmax; 933 indicies total

         0 seconds required to construct component=gen_pmin; 933 indicies total

         0 seconds required to construct component=gen_heatrate; 933 indicies total

         0 seconds required to construct component=gen_fuelcost; 933 indicies total

         0 seconds required to construct component=gen_sucost; 933 indicies total

         0 seconds required to construct component=gen_sdcost; 933 indicies total

         0 seconds required to construct component=gen_nlcost; 933 indicies total

         0 seconds required to construct component=gen_mindowntime; 933 indicies total

         0 seconds required to construct component=gen_minuptime; 933 indicies total

         0 seconds required to construct component=gen_hourlyramp; 933 indicies total

         0 seconds required to construct component=gen_spramp; 933 indicies total

         0 seconds required to construct component=gen_nspramp; 933 indicies total

         0 seconds required to construct component=gen_susdramp; 933 indicies total

         0 seconds required to construct component=gen_zone1; 933 indicies total

         0 seconds required to construct component=gen_zone2; 933 indicies total

         0 seconds required to construct component=gen_zone3; 933 indicies total

         0 seconds required to construct component=gen_zone4; 933 indicies total

         0 seconds required to construct component=res_req; 24 indicies total

         0 seconds required to construct component=res_req_zone1; 24 indicies total

         0 seconds required to construct component=res_req_zone2; 24 indicies total

         0 seconds required to construct component=res_req_zone3; 24 indicies total

         0 seconds required to construct component=res_req_zone4; 24 indicies total

      0.29 seconds required to construct component=hourly_load; 103056 indicies total

         0 seconds required to construct component=hourly_wind; 1848 indicies total

      0.06 seconds required to construct component=gen_status; 22392 indicies total

      0.06 seconds required to construct component=gen_vstatus; 22392 indicies total

      0.06 seconds required to construct component=gen_wstatus; 22392 indicies total

      0.08 seconds required to construct component=gen_contingency; 22392 indicies total

      0.86 seconds required to construct component=branch_contingency; 216000 indicies total

         0 seconds required to construct component=net_interchange; 24 indicies total

         0 seconds required to construct component=demand_multiplier; 24 indicies total

         0 seconds required to construct component=vang; 1 indicies total

    235.30 seconds required to construct component=ptdf; 69534000 indicies total

         0 seconds required to construct component=masterU; 0 indicies total

         0 seconds required to construct component=masterV; 0 indicies total

         0 seconds required to construct component=masterW; 0 indicies total

      0.11 seconds required to construct component=gen_supply; 22392 indicies total

      0.11 seconds required to construct component=sp_res; 22392 indicies total

      0.11 seconds required to construct component=nsp_res; 22392 indicies total

      0.23 seconds required to construct component=injection; 185424 indicies total

      1.11 seconds required to construct component=negativeviolation; 216000 indicies total

      1.11 seconds required to construct component=positiveviolation; 216000 indicies total

         0 seconds required to construct component=c_index; 0 indicies total

         0 seconds required to construct component=c; 0 indicies total

      3.91 seconds required to construct component=obj; 1 indicies total

             Cloning detected! (clone counters: 2, 0, 0)

      556.05 seconds required to construct instance=unknown

        1.23 seconds required for preprocessing

The model stalls here for a very long time before constructing the associated constraints. 

I have tried modelling all the constraints as a constraint list as well as without using constraint list. Thank you.

Regards,

Pranav

[Gabriel Hackebeil]

Is this the large instance or the small instance? And does is “stall for a very long time” before you see this output or after?

Gabe

[pranav]

Hi Gabe,

Thank you for the suggestions. I am using a .dat file. But I don't think the problem is with reading the data file. I guess the problem is with constructing the constraint rules.

Regards,

Pranav

On Thursday, August 6, 2015 at 7:18:02 AM UTC-6, pranav wrote:

[pranav]

This is for large instance. It stalls for very long time after I see this output. I have modeled the constraints only after this stage and it does not build any constraints.

Regards,

Pranav

[pranav]

Hi Cristian,

Yes I guess the problem is with the construction of the constraints. Whether I specify the constraints as rules or as constraint list using  'for' lops, it does not help. Thank you.

Regards,

Pranav 

[Gabriel Hackebeil]

This is for large instance. It stalls for very long time after I see this output. I have modeled the constraints only after this stage and it does not build any constraints.

Okay, that makes more sense.

Yes I guess the problem is with the construction of the constraints. Whether I specify the constraints as rules or as constraint list using  'for' lops, it does not help. Thank you.

I think what Cristian means is that the way you construct the constraint expressions can significantly impact construction time. For instance:

expr = 0.0

for i in model.some_set:

     expr = expr + model.x[i]

return expr

Is much much slower than

return sum(model.x[i] for i in model.some_set[i])

Also, as I said before, PyPy is definitely worth looking into. You can typically execute identical Python code an order of magnitude faster than with the standard CPython implementation.

Gabe

[Watson, Jean-Paul]

There are a few things we can try, but let’s look at the 73-bus system first, specifically in terms of instantiation time. Assuming you are (or can) run the “pyomo” script with the additional option “—report-timing”, this will display information as follows:

[    0.00] Setting up Pyomo environment

[    0.00] Applying Pyomo preprocessing actions

[    0.00] Creating model

         0 seconds required to construct component=FOOD; 9 indicies total

         0 seconds required to construct component=cost; 9 indicies total

         0 seconds required to construct component=f_min; 9 indicies total

         0 seconds required to construct component=f_max; 9 indicies total

         0 seconds required to construct component=NUTR; 7 indicies total

         0 seconds required to construct component=n_min; 7 indicies total

         0 seconds required to construct component=n_max; 7 indicies total

         0 seconds required to construct component=amt_index; 63 indicies total

         0 seconds required to construct component=amt; 63 indicies total

         0 seconds required to construct component=Buy; 9 indicies total

         0 seconds required to construct component=Total_Cost; 1 indicies total

         0 seconds required to construct component=Entree; 1 indicies total

         0 seconds required to construct component=Side; 1 indicies total

         0 seconds required to construct component=Drink; 1 indicies total

        0.00 seconds required to construct instance=unknown

        0.00 seconds required for problem transformations

[    0.10] Applying solver

[    0.13] Processing results

    Number of solutions: 1

    Solution Information

      Gap: 0.0

      Status: optimal

      Function Value: 2.81

Of course, this isn’t particularly interesting for the diet problem, but it should help to identify which components in your model are taking the most time. 

We would also be interested in running this locally; if you are willing to send to “jwa...@sandia.gov”, I can circulate so that we can have a look.

jpw

From: <pyomo...@googlegroups.com> on behalf of pranav <pranava...@gmail.com>
Reply-To: "pyomo...@googlegroups.com" <pyomo...@googlegroups.com>
Date: Thursday, August 6, 2015 at 7:18 AM
To: "pyomo...@googlegroups.com" <pyomo...@googlegroups.com>
Subject: [EXTERNAL] Excessively long time to create models for large systems

--

[Watson, Jean-Paul]

Sorry for answering a bit out of order, but…

I suspect the problem is that you are constructing a dense branch-by-bus incidence matrix – this is order (m times n), which his large. In reality, this is quite sparse, as I’m sure you know. So, what you really want to have is a set of pairs that specify the indices that are valid/relevant. 

Does this make sense?

There are others on this list (Dave, in particular) that might be able to dig up and send you a nice set of slides on how to model such sparse sets and associated constraints.

jpw

From: <pyomo...@googlegroups.com> on behalf of Gabriel Hackebeil <gabe.ha...@gmail.com>
Reply-To: "pyomo...@googlegroups.com" <pyomo...@googlegroups.com>
Date: Thursday, August 6, 2015 at 7:56 AM
To: "pyomo...@googlegroups.com" <pyomo...@googlegroups.com>

[pranav]

Thank you JPW. I would like to go over the slides for modeling sparse sets. It would be really helpful if you can share it with me.

Regards,

Pranav

[pranav]

Hi JPW, 

The constraints corresponding to the dense branch by bus matrix took the maximum time in the 73 bus system. However, it took less than 2 seconds. In the larger instance case it took 3 hours to construct 1 set of constraint corresponding to the branch by bus matrix.

Regards,

Pranav

On Thursday, August 6, 2015 at 11:49:06 AM UTC-6, JPW wrote:

[David Woodruff]

Really, the slides are just the slides version of the text in the section called Sparse Index Sets in the documentation at pyomo.org

  Dave

[Watson, Jean-Paul]

Pranav: Please take a look at the slides Dave references, and let us know if you have any questions. Sparsifying the bus/branch incidence matrix should mitigate most of the run time issues you are observing.

jpw

[Watson, Jean-Paul]

I can’t tell from your model the specific typeof BRANCHxBUS, but I’m assuming it is a set, no? That’s taking 311 seconds to construct, presumably because it is dense. I’m assuming (?) that ptdf is a parameter? If so, it’s also taking a huge chunk of time due to the density – there is no way you have 69M of anything in a 76 bus system! (?).

jpw

From: <pyomo...@googlegroups.com> on behalf of pranav <pranava...@gmail.com>
Reply-To: "pyomo...@googlegroups.com" <pyomo...@googlegroups.com>

Date: Thursday, August 6, 2015 at 1:17 PM
To: "pyomo...@googlegroups.com" <pyomo...@googlegroups.com>

[David Woodruff]

I have attached the slides. They are a little old (e.g., they are from the coopr days) so I recommend looking at the Sparse Index Sets section of the online documentation at Pyomo.org

[pranav]

Hi Dave,

Thanks for the reference. 

JPW, 

The ptdf is given as a parameter. The branch and bus are sets. BRANCHxBUS is constructed as given below:

def branchbusrule(model):

ans = set()

for k in model.branch:

for n  in model.bus:

ans.add((k,n))

return ans

model.BRANCHxBUS = Set(initialize=branchbusrule, dimen=2)

All the information I have sent out earlier are for 9000 bus model. Its not for 73 bus system.

Regards,

Pranav

[Watson, Jean-Paul]

On the PTDF front, I think you want a filter in the doubly nested loop that only adds the tuple if the bus and branch are actually connected. You’ll still end up doing a lot of unnecessary work, and things would go quicker still if you had something that defined adjacency in your .dat file (or somewhere in data, if you are using concrete models). 

So the timing results are for the full 9K bus system? This seems to be running in about 500 seconds, in terms of model instantiation. With the above modifications, that time will go to nil. Am I missing something?

[Pranavamoorthy Balasubramanian]

Hi JPW,

It has not constructed the constraints yet. For constructing a constraint like the one given below, it takes close to 3 hours and close to 26 GB RAM memory. I have given all the constraints as constraint list. the constraint list is created after creating the instance.

instance = model.create(datafile+ '.dat', report_timing=True)

for (k,t) in instance.BRANCHxTIME:

if instance.branch_contingency[k,t]!=1:

instance.c.add(sum(instance.ptdf[k,n]*instance.injection[n,t] for n in instance.bus) + instance.negativeviolation[k,t] >= -instance.branch_maxcap[k])

After creating all constraints, I do preprocess()

opt = SolverFactory("cbc")

opt.options['threads'] = 16

results = opt.solve(instance, suffixes=['rc'], tee=True)

instance.load(results)

Thank you.

Regards,

Pranav

--
You received this message because you are subscribed to a topic in the Google Groups "Pyomo Forum" group.
To unsubscribe from this topic, visit https://groups.google.com/d/topic/pyomo-forum/YTPByhZzboA/unsubscribe.
To unsubscribe from this group and all its topics, send an email to pyomo-forum...@googlegroups.com.

For more options, visit https://groups.google.com/d/optout.

--

Thanks and Regards

[Gabriel Hackebeil]

I just want to chime in here and say that I had trouble figuring out how to define a sparse adjacency set through DAT files when I first started. So here are some tips:

(1) We are saying that there are likely more indices of ptdf that are zero than there are not, so don’t even define a parameter or an index set that is BRANCH x BUS in size, only define a parameter and index set over the nonzero elements of ptdf (and do this in whatever code is creating the data files). If you don’t believe us, please report the number of indices of ptdf that are nonzero (out of the total 69534000 indices).

(2) You are going to define a Set on your model that looks something like “model.BRANCHxBUS_sparse = Set(dimen=2)”. Your DAT file is going to look something like:

set BRANCH := b1 b2 b3 b5 b6 b7;

set BUS := s1 s2;

set BRANCHxBUS_sparse := (b1,s1) (b2,s1) (b2,s2) (b3,s2) (b7,s1);

where BRANCHxBUS_sparse could get pretty long if the sparse set is large (but it should have much much fewer than |BRANCH x BUS| elements).

(3) If you want to skip the DAT file format and work with a ConcreteModel, you can try storing the data for ptdf in Matrix Market format (google that) and then use the Python Scipy package to load the sparse matrix into memory (http://docs.scipy.org/doc/scipy/reference/generated/scipy.io.mmread.html#scipy.io.mmread). This will likely lead to faster load times if your sparse set is large. It’s probably the better route if ptdf is dense as well (which hopefully it is not).

Gabe

[Pranavamoorthy Balasubramanian]

Thank you for the information Gabe. I will try it out.

Regards,

Pranav

[Watson, Jean-Paul]

Let’s first try to get the sets defined as Gabe suggested, and then we’ll move to worrying about the constraints. 

jpw

[Amandeep Gupta]

Hi,

    I am trying to build an economic dispatch program for 118 bus system using Pyomo. I have been facing similar problems as Pranav did with model build-up times. I have followed the suggestions on this thread to reduce the time for component construction. I am, however, stuck at the constraint building stage. Pyomo is taking a lot of time to build line flow constraints with the ptdf matrix. I have attached the report-timing output, I realize that I am doing something wrong, hence, the cloning detected warning. I am not sure how to deal with it though. Can you please give some suggestions?

Thanks

Aman

[    0.00] Setting up Pyomo environment

[    0.00] Applying Pyomo preprocessing actions

[    0.01] Creating model

         0 seconds required to construct component=GEN; 54 indicies total

         0 seconds required to construct component=nodes; 118 indicies total

         0 seconds required to construct component=nodes_volt; 117 indicies total

         0 seconds required to construct component=lines; 179 indicies total

         0 seconds required to construct component=horizon; 1 indicies total

         0 seconds required to construct component=hours; 24 indicies total

         0 seconds required to construct component=twoplushours; 23 indicies total

         0 seconds required to construct component=threeplushours; 22 indicies total

         0 seconds required to construct component=variab_cost; 54 indicies total

         0 seconds required to construct component=RU_index; 1296 indicies total

      0.01 seconds required to construct component=RU; 1296 indicies total

         0 seconds required to construct component=RD_index; 1296 indicies total

      0.01 seconds required to construct component=RD; 1296 indicies total

         0 seconds required to construct component=maxproduce_lim; 54 indicies total

         0 seconds required to construct component=minproduce_lim; 54 indicies total

         0 seconds required to construct component=line_cap; 179 indicies total

         0 seconds required to construct component=bus_gen_matr_index; 6372 indicies total

      0.02 seconds required to construct component=bus_gen_matr; 6372 indicies total

         0 seconds required to construct component=bus_line_matr_index; 21122 indicies total

      0.07 seconds required to construct component=bus_line_matr; 21122 indicies total

         0 seconds required to construct component=bl_index; 21122 indicies total

      0.07 seconds required to construct component=bl; 21122 indicies total

         0 seconds required to construct component=demand_index; 2832 indicies total

      0.01 seconds required to construct component=demand; 2832 indicies total

         0 seconds required to construct component=unct_nodes; 0 indicies total

         0 seconds required to construct component=diagonal_D_index; 32041 indicies total

      0.10 seconds required to construct component=diagonal_D; 32041 indicies total

         0 seconds required to construct component=shift_factor_index; 21122 indicies total

      0.06 seconds required to construct component=shift_factor; 21122 indicies total

         0 seconds required to construct component=produce_p_index; 1296 indicies total

         0 seconds required to construct component=produce_p; 1296 indicies total

         0 seconds required to construct component=line_trans_index; 4296 indicies total

      0.01 seconds required to construct component=line_trans; 4296 indicies total

         0 seconds required to construct component=theta_index; 2808 indicies total

         0 seconds required to construct component=theta; 2808 indicies total

         0 seconds required to construct component=pinject_index; 2832 indicies total

         0 seconds required to construct component=pinject; 2832 indicies total

         0 seconds required to construct component=genstatus_x_index; 1296 indicies total

         0 seconds required to construct component=genstatus_x; 1296 indicies total

         0 seconds required to construct component=slack_lines_index; 4296 indicies total

      0.01 seconds required to construct component=slack_lines; 4296 indicies total

         0 seconds required to construct component=slack_demand1; 24 indicies total

         0 seconds required to construct component=slack_demand2; 24 indicies total

         0 seconds required to construct component=rampup_index; 1296 indicies total

      0.02 seconds required to construct component=rampup; 1296 indicies total

         0 seconds required to construct component=rampdown_index; 1296 indicies total

      0.03 seconds required to construct component=rampdown; 1296 indicies total

         0 seconds required to construct component=flow_upper_bound_index; 4296 indicies total

      0.07 seconds required to construct component=flow_upper_bound; 4296 indicies total

         0 seconds required to construct component=flow_lower_bound_index; 4296 indicies total

      0.07 seconds required to construct component=flow_lower_bound; 4296 indicies total

         0 seconds required to construct component=gen_upper_bound_index; 1296 indicies total

      0.02 seconds required to construct component=gen_upper_bound; 1296 indicies total

         0 seconds required to construct component=gen_lower_bound_index; 1296 indicies total

      0.02 seconds required to construct component=gen_lower_bound; 1296 indicies total

      0.01 seconds required to construct component=energy_conserv; 24 indicies total

         0 seconds required to construct component=flow_upper_bound1_index; 4296 indicies total

    220.20 seconds required to construct component=flow_upper_bound1; 4296 indicies total

             Cloning detected! (clone counters: 217680)

      0.07 seconds required to construct component=opt_mincost_mstr; 1 indicies total

             Cloning detected! (clone counters: 218976)

      220.88 seconds required to construct instance=unknown

        0.00 seconds required for problem transformations

[  222.67] Applying solver

        2.13 seconds required to write file

        2.14 seconds required for presolve

        0.37 seconds required for solver

        0.00 seconds required to read logfile 

        0.20 seconds required to read solution file 

        0.30 seconds required for postsolve

[  225.48] Processing results

    Number of solutions: 1

    Solution Information

      Gap: 0.0

      Status: optimal

      Function Value: 2272732.93679

    Solver results file: results.json

[  225.78] Applying Pyomo postprocessing actions

[  225.78] Pyomo Finished

[Watson, Jean-Paul]

Hi Aman,

 

Would you be willing to share your model, off-forum? Or at least the code associated with flow_upper_bound1? Cloning is the likely culprit, and it should be able to eliminate that by re-working implementation specifics of that rule.

 

jpw

[Bill Hart]

I would definitely try using Python 3.6.  There are known performance issues in Python 3.4 that might impact your analysis.

--Bill

On Thursday, August 6, 2015 at 7:18:02 AM UTC-6, pranav wrote:

[Bill Hart]

@Gabe:

The current Pyomo release does not support PyPy robustly.  We won't have a Pyomo release supporting PyPy until later this spring.

--Bill

To unsubscribe from this group and stop receiving emails from it, send an email to pyomo-forum+unsubscribe@googlegroups.com.

[Bill Hart]

@JP, @Aman:

I'd also like a copy of the model.  My rework of the expression system has eliminated the excessive cloning that can be done in the latest Pyomo release.  But I'd like to confirm this.

--Bill

To unsubscribe from this group and stop receiving emails from it, send an email to pyomo-forum+unsubscribe@googlegroups.com.
For more options, visit https://groups.google.com/d/optout.

[Watson, Jean-Paul]

I’ll leave it to @Aman as to whether he’s willing to share the model / data. In our off-line discussions, he managed to (largely) mitigate the model creation bottleneck by using sparse representations of what were being represented as dense matrices (with a lot of 0s). So while it would be a good test case for profiling purposes…

 

jpw

[Siirola, John D]

I should also point out that Pyomo 4.0.9638 is 3 years old.  There have been numerous improvements over that time.  Please upgrade to the current release (Pyomo 5.3) before reporting performance problems.

 

In particular there is *no* support for PyPy in 4.0.  You should be on a 5.x release before even *trying* pypy. 

 

john

 

From: pyomo...@googlegroups.com [mailto:pyomo...@googlegroups.com] On Behalf Of Bill Hart
Sent: Wednesday, February 07, 2018 12:54 PM
To: Pyomo Forum <pyomo...@googlegroups.com>

To unsubscribe from this group and stop receiving emails from it, send an email to pyomo-forum...@googlegroups.com.

For more options, visit https://groups.google.com/d/optout.

 

--

You received this message because you are subscribed to the Google Groups "Pyomo Forum" group.

To unsubscribe from this group and stop receiving emails from it, send an email to pyomo-forum...@googlegroups.com.

[serg]

Hi,
I am also having performance issues with medium-scale models. I have identified the two problematic constraints using profiling. They both involve highly sparse matrices so I do not know if that could be the issue How can I go about optimizing these constraints? I could not find any documentation online, other than this thread.
-Sergio.

[Watson, Jean-Paul]

Performance can be significantly impacted by the way in which you specify the matrices – constraint coefficients, in particular. If you use what is effectively a dense representation, e.g., with Param(x,y,default=0.0), then you will see very long instantiation times.

 

If you share the constraint rules and any definitions of associated sets and parameters, we can say more.

 

jpw

[serg]

Thank you for your response. I am really looking forward to using Pyomo, but currently the performance is a limitation, I am most likely not coding something correctly. I looked into sparse sets, but I was not able to apply it for the sparse parameter matrix S. Here are the problematic constraints:

def mass_balance(model, i, k, m):
    if not any(model.S[i, j, k] for j in model.J):
        return Constraint.Skip
    return sum(model.S[i, j, k] * model.v[j, k, m] for j in model.J) == model.b[i, k]
model.mass_balance = Constraint(model.I, model.K, model.M, rule=mass_balance)

def dual_con(model, j, k, m):
    return sum(model.lambda1[i, k, m] * model.S[i, j, k] for i in model.I) \
           + model.muu[j, k, m] - model.mul[j, k, m] == model.c[j, k, m]
 
model.dual_con = Constraint(model.J, model.K, model.M, rule=dual_con)

With variable definitions:
model.I = Set()
model.J = Set()
model.K = Set()
model.M = Set()

model.S = Param(model.I, model.J, model.K, default=0)

To unsubscribe from this group and stop receiving emails from it, send an email to pyomo-forum+unsubscribe@googlegroups.com.
For more options, visit https://groups.google.com/d/optout.

[Siirola, John D]

The problem with your model is that while model.S is stored sparsely, it is technically a dense parameter (there are valid values for all indices).  Your constraint is actually working over the dense matrix (twice!).

 

Something like the following would be faster (but still dense):

 

def mass_balance(model, i, k, m):

    flows = sum(model.S[i, j, k] * model.v[j, k, m] for j in model.J if model.S[i,j,k] != 0)

    if flows is 0 or flows is 0.0:

        return Constraint.Skip

    return (flows, model.b[i, k])  # you could do return flows == model.b[I,k], but it would be slightly slower in the current system

model.mass_balance = Constraint(model.I, model.K, model.M, rule=mass_balance)

def dual_con(model, j, k, m):

    return sum(model.lambda1[i, k, m] * model.S[i, j, k] for i in model.I if model.S[i,j,k] != 0) \

           + model.muu[j, k, m] - model.mul[j, k, m] == model.c[j, k, m]
model.dual_con = Constraint(model.J, model.K, model.M, rule=dual_con)

The “right” solution is to use a data structure that directly exposes the relevant sparsity.  If you know the sparsity, great!  If you need to infer it from the data, then something like this would work (putting it in the model after model.S is declared, of course):

 

model.SPARSE_IK = Set(model.I, model.K)

model.SPARSE_JK = Set(model.J, model.K)

 

def detect_sparsity(model):

    for i,j,k in model.S.sparse_iterkeys():

        model.SPARSE_IK[i,k].add(j)

        model.SPARSE_JK[j,k].add(i)

model.detect_sparsity = BuildAction(rule=detect_sparsity)

 

def mass_balance(model, i, k, m):

    if len(model.SPARSE_IK[i,k]) == 0:

        return Constraint.Skip

    return sum(model.S[i, j, k] * model.v[j, k, m] for j in model.SPARSE_IK[i,k]) == model.b[i,k]

model.mass_balance = Constraint(model.I, model.K, model.M, rule=mass_balance)

def dual_con(model, j, k, m):

    return sum(model.lambda1[i, k, m] * model.S[i, j, k] for i in model.SPARSE_JK[j,k]) \

           + model.muu[j, k, m] - model.mul[j, k, m] == model.c[j, k, m]
model.dual_con = Constraint(model.J, model.K, model.M, rule=dual_con)

The “BuildAction” is so that this will work for both Abstract and Concrete models.  If you are in Concrete land, you can avoid it by just running (or inlining) the function.

 

Assuming I didn’t make any typos, this should scale with the number of nonzeros (and not the size of the indices).

 

john

[serg]

Thank you! It works! Now I can actually build a large model in ~22 seconds, and the current bottleneck is writing to .lp which takes ~40s. That is a minute though, so it is completely acceptable. Thanks again.

To unsubscribe from this group and stop receiving emails from it, send an email to pyomo-forum+unsubscribe@googlegroups.com.

For more options, visit https://groups.google.com/d/optout.

[Michael]

Hi all,

I'm relatively new to pyomo and I'm also having a similar issue. I have tried to use the advice from above but I am struggling to implement it, (I have constructed the for loops using the method above, but even after reading the documentation I'm not sure about sparse sets).

I have some constraints that are taking a long time to construct (about 1min just for one constraint, I have many constraints so it takes a while to solve the whole model). My model is concrete and I am using dictionaries to initialize the params. I have attached an example code below to help explain my task. The example constraint is only medium-sized, I have others that are about twice the size. I would be very grateful for any help or suggestions that can help speed up my model.

[Dennis Gray]

I have a feeling that it’s just slow to construct your model. For example, I’m building a spatial model in pyomo and experimenting with different grid sizes. When I jumped from making a 20x20 grid for 10 time periods to a 30x30 grid for 10 time periods, I go from making 4000 constraints to 9000 constraints and it significantly increased model creation time. 

I would suggest just letting pyomo do its thing, but before you initiate the solving command (if you’re scripting), you can input some code to indicate when all of the constraints have been created by something like a simple statement of your choice. 

--
You received this message because you are subscribed to the Google Groups "Pyomo Forum" group.
To unsubscribe from this group and stop receiving emails from it, send an email to pyomo-forum...@googlegroups.com.

To view this discussion on the web visit https://groups.google.com/d/msgid/pyomo-forum/adc3c3a6-79d3-4f27-a6d8-0817938cf88f%40googlegroups.com.

--

Cheers,

Dennis Gray

MSc Candidate - Department of Environmental Economics and Resource Sociology

116 st and 85 ave, Edmonton AB Canada, T6G 2R3

C: 780-717-4190