Recent Questions - Operations Research Stack Exchange most recent 30 from or.stackexchange.com 2021-11-29T03:23:55Z https://or.stackexchange.com/feeds https://creativecommons.org/licenses/by-sa/4.0/rdf https://or.stackexchange.com/q/7348 3 Range limits on terms in the objective function of an LP Henry https://or.stackexchange.com/users/2564 2021-11-27T16:50:05Z 2021-11-28T12:22:59Z <p>I have a linear maximization problem with an objective as follows: <span class="math-container">$$\sum c_i\cdot\text{exp}_i$$</span> where <span class="math-container">$c_i$</span> are constants (positive or negative) and <span class="math-container">$\text{exp}_i$</span> are linear expressions of the free variables, which may be positive or negative.</p> <p>I want to put a floor and ceiling on the <code>exp[i]</code> in the objective, so that no term may add or subtract more than <code>abs(c[i] * threshold)</code> to/from the objective value.</p> <p>It's straight-forward to implement the ceiling for <code>c[i] &gt; 0</code> and the floor for <code>c[i] &lt; 0</code> as penalties:</p> <pre><code>if c[i] &gt; 0: penalty[i] &gt; exp[i] - threshold penalty[i] &gt; 0 subtract c[i] * penalty[i] from the objective function if c[i] &lt; 0: penalty[i] &gt; threshold - exp[i] penalty[i] &gt; 0 subtract c[i] * penalty[i] from the objective function </code></pre> <p>However this method does not work for a floor when <code>c[i] &gt; 0</code> or for a ceiling when <code>c[i] &lt; 0</code> because <code>penalty[i]</code> is not bounded from above.</p> <p>So I am asking if there is a different &quot;trick&quot; that can be applied here, or can this even be done as a linear problem, or alternately is there heuristic approach etc.</p> https://or.stackexchange.com/q/7347 0 How to solve this linear programming minimization problem with the BigM method TheInquirer https://or.stackexchange.com/users/7612 2021-11-26T19:25:21Z 2021-11-27T08:24:20Z <p>Given this minimization problem: <span class="math-container">$$z=100x_1+77x_2+80x_3\\10x_1+7x_2+2x_3 \geq 12\\2x_1+3x_2+4x_3 \geq 3\\x_1+2x_2+x_3 \geq 1\\x_1,x_2,x_3 \geq 0$$</span> After adding artificial variables: <span class="math-container">$$z-100x_1-77x_2-80x_3-MR_1-MR_2-MR_3=0\\10x_1+7x_2+2x_3-s_1+R_1=12\\2x_1+3x_2+4x_3-s_2+R_2=3\\x_1+2x_2+x_3-s_3+R_3=1$$</span> What do you do from here?</p> https://or.stackexchange.com/q/7343 9 Beginner friendly open source projects in O.R PSLP https://or.stackexchange.com/users/405 2021-11-26T14:51:16Z 2021-11-26T18:09:25Z <p>The title says it all: <strong>Are there any open source projects, that are open to new contributers?</strong></p> <ul> <li>If not, should there be?</li> <li>If yes, what are they and what is the best way to get involved?</li> </ul> https://or.stackexchange.com/q/7338 1 Check if there exists a vector that satisfies a set of inequalities but violates another set of inequalities Duh Huh https://or.stackexchange.com/users/6094 2021-11-24T17:40:49Z 2021-11-24T19:40:11Z <p><strong>Problem</strong></p> <p>Given rectangular matrices <span class="math-container">$A$</span>, <span class="math-container">$B$</span> and vectors <span class="math-container">$\vec{a}$</span>, <span class="math-container">$\vec{b}$</span>, how to check if there exists an <span class="math-container">$\vec{x}$</span> that satisfies the following conditions?</p> <ol> <li><span class="math-container">$\vec{x} \succcurlyeq \vec{0}$</span> is true.</li> <li><span class="math-container">$A\vec{x} \succcurlyeq \vec{a}$</span> is true.</li> <li><span class="math-container">$B\vec{x} \succcurlyeq \vec{b}$</span> is false.</li> </ol> <p>For vectors <span class="math-container">$\vec{u}$</span> and <span class="math-container">$\vec{v}$</span>, <span class="math-container">$\vec{u} \succcurlyeq \vec{v}$</span> means:</p> <ul> <li>The vectors have the same number of elements.</li> <li>For all i, the i<sup>th</sup> elements of the vectors satisfy <span class="math-container">$u_i &gt; v_i$</span>.</li> </ul> <p>I can't use <span class="math-container">$B\vec{x} \preccurlyeq \vec{b}$</span> because <span class="math-container">$\vec{x}$</span> only have to violate at least 1 inequalities among all inequalities represented by <span class="math-container">$B\vec{x} \succcurlyeq \vec{b}$</span>.</p> <p>Each problem has about 100 variables and 150 inequalities. I want to solve at least 10,000 problems per minute on a laptop. All of the problems share a large number of inequalities.</p> <hr /> <p><strong>Ideas</strong></p> <p>Define two problems:</p> <ul> <li><span class="math-container">$P_{A}$</span> means <span class="math-container">$\vec{x} \succcurlyeq \vec{0}$</span> and <span class="math-container">$A\vec{x} \succcurlyeq \vec{a}$</span> are true.</li> <li><span class="math-container">$P_{B}$</span> means <span class="math-container">$\vec{x} \succcurlyeq \vec{0}$</span> and <span class="math-container">$B\vec{x} \succcurlyeq \vec{b}$</span> are true.</li> </ul> <p>If I remember correctly, I think:</p> <ul> <li>Vectors that satisfies <span class="math-container">$P_{A}$</span> are the convex combinations of the vertex set <span class="math-container">$V_{A}$</span>. <ul> <li>The vectors in <span class="math-container">$V_A$</span> form a matrix <span class="math-container">$M_{A}$</span>.</li> <li>The convex combinations form a convex set <span class="math-container">$S_{A}$</span>.</li> </ul> </li> <li>Vectors that satisfies <span class="math-container">$P_{B}$</span> are the convex combinations of the vertex set <span class="math-container">$V_{B}$</span>. <ul> <li>The vectors in <span class="math-container">$V_B$</span> form a matrix <span class="math-container">$M_{B}$</span>.</li> <li>The convex combinations form a convex set <span class="math-container">$S_{B}$</span>.</li> </ul> </li> <li>Intersection of two convex sets is convex.</li> <li>Difference of two convex sets is not necessarily convex.</li> </ul> <p><strong>Idea 2</strong></p> <ul> <li>The problem is equivalent to showing if <span class="math-container">$S_{A}$</span> is not a subset of <span class="math-container">$S_{B}$</span>.</li> <li>I think the problem is also equivalent to showing if there exists a vertex <span class="math-container">$\vec{\alpha} \in V_A$</span> such that <span class="math-container">$\vec{\alpha}$</span> is not a convex combination of <span class="math-container">$V_B$</span>. That means <span class="math-container">$M_{B}\vec{u} = \vec{\alpha}$</span> either has no solution or <span class="math-container">$\vec{u}\cdot\vec{1} \ne 1$</span>.</li> </ul> <p>I don't know how to find the vertexes. And there are potentially so many vertexes.</p> <p><strong>Idea 3</strong></p> <ul> <li>Start from an old problem without such a <span class="math-container">$\vec{x}$</span></li> <li>By comparing the new <span class="math-container">$P_{B}$</span> and old <span class="math-container">$P_{B}$</span>, we can identify the new cuts.</li> <li>Invert the comparison in the new cuts. (<span class="math-container">$\ge$</span> to <span class="math-container">$\le$</span>).</li> <li>Check if each new inverted cut is compatible with the new <span class="math-container">$P_{A}$</span>. <ul> <li>For each inverted new cut, add the inverted new cut to the new <span class="math-container">$P_A$</span>. If simplex algorithm says the combined problem is feasible, the new inverted cut is compatible with the new <span class="math-container">$P_A$</span>.</li> </ul> </li> <li>The new problem has such a <span class="math-container">$\vec{x}$</span> if and only if we find a new inverted cut that is incompatible with the new <span class="math-container">$P_{A}$</span>.</li> </ul> <p>That would be at most 60 runs of simplex per problem. The method caches the inequalities instead of the vertices.</p> <p>I don't need the solution and only need to know if the new inverted cut is compatible with the existing cuts. Is there a faster way to get this?</p> https://or.stackexchange.com/q/7336 3 If variable falls below a certain value, include difference to set value in objective Balasar https://or.stackexchange.com/users/7594 2021-11-24T15:03:44Z 2021-11-24T17:16:11Z <p>I think its easiest to describe my goal first and continue with my implementation and the resulting problems!</p> <p>My goal: Using Pyomo as interface and Gurobi as solver, if a variable <span class="math-container">$x_{i,t}$</span> falls below a certain value <span class="math-container">$b_{i,t}$</span> the deviation to given threshold shall be added to the objective, but only if it is below. <span class="math-container">$i$</span> denotes a state inside a state space representation and <span class="math-container">$t$</span> the current time step. So my objective should describe something like <span class="math-container">$$J=\sum_i^S\sum_t^T(\mathrm{if}\, x_{i,t}&lt;b_{i,t}:b_{i,t}-x_{i,t}, \mathrm{else}\,0).$$</span> So it is comparable to a ReLU activation function</p> <p>My Implementation/Problem: I tried implementing a binary variable <span class="math-container">$y_{i,t}\in\{0;1\}$</span> using BigM method <span class="math-container">$$b_{i,t}-x_{i,t}\le My_{i,t}$$</span> and describing <span class="math-container">$$J=\sum_i^S\sum_t^Ty_{i,t}\cdot(b_{i,t} - x_{i,j}).$$</span> The problem is if <span class="math-container">$x_{i,t}$</span> is above the threshold the optimizer also sets <span class="math-container">$y_{i,t}=1$</span> since this allows for negative objective terms, which minimizes the objective but is not what I intended.</p> <p>Is there any way on how to implement this? In addition, later I want to implement the same just with an upper threshold. Both of these thresholds should be usable independently of each other (only the upper or only the lower), but also in combination with each other, creating a band where no objective is added (this is probably just an summation of both onbjectives).</p> <p>Maybe as an additional note, I should mention that <span class="math-container">$x_{i,t+1}$</span> is calculated using Pyomo DAE based on its differential <span class="math-container">$\dot{x}_{i,t}$</span>, the systems input <span class="math-container">$u_{j,t}$</span> and the internal integration scheme of Pyomo DAE (so basically a state stace). The first value <span class="math-container">$x_{i,0}$</span> is fixed via an equality constraint so the optimizer can only really manipulate <span class="math-container">$u$</span> since the differential at a given time step <span class="math-container">$\dot{x}_{i,t}$</span> is calculated via the current values of <span class="math-container">$x_{i,t}$</span> and <span class="math-container">$u_{j,t}$</span>. Both <span class="math-container">$x_{i,t}$</span> and <span class="math-container">$u_{j,t}$</span> are bounded but not necessarily <span class="math-container">$\ge 0$</span>.</p> https://or.stackexchange.com/q/7334 6 Network design problem with rounded capacity constraints Joris Kinable https://or.stackexchange.com/users/49 2021-11-24T02:13:03Z 2021-11-24T08:51:55Z <p>I have a network design problem with complicating capacity constraints which I'm trying to model through a Mixed Integer Programming formulation.</p> <p>The problem is defined on a directed, incomplete graph <span class="math-container">$G(V,A)$</span>. A binary variable <span class="math-container">$x_{uv}^k$</span> defines whether commodity <span class="math-container">$k\in K$</span> is routed via arc <span class="math-container">$(u,v)\in A$</span>. Parameter <span class="math-container">$q_k\in \mathbb{R}_{&gt;0}$</span> defines the volume of commodity <span class="math-container">$k\in K$</span>. For a subset of nodes <span class="math-container">$U\subset V$</span>, I have the following capacity constraints:</p> <p><span class="math-container">\begin{align} &amp; \sum_{v:(u,v)\in A}\Big\lceil \sum_{k\in K}q_kx_{uv}^k\Big\rceil_{\ell}\leq Q_u &amp; \forall u\in U \end{align}</span> Here, <span class="math-container">$Q_u$</span> is the capacity of node <span class="math-container">$u\in U$</span>, and <span class="math-container">$\lceil\cdot\rceil_{\ell}$</span> is a rounding operator that rounds <strong>up</strong> to the nearest multiple of <span class="math-container">$\ell$</span>. In my application, <span class="math-container">$\ell$</span> can take the values <span class="math-container">$0.5$</span> or <span class="math-container">$1$</span>.</p> <p>The rounding operation makes it troublesome to formulate this constraint. In order to model this, I could associate non-negative <strong>integer</strong> helper variables <span class="math-container">$p_{uv}\in \mathbb{Z}_{\geq 0}$</span> with all arcs <span class="math-container">$(u,v)\in A$</span>, and then state the following two constraints: <span class="math-container">\begin{align} &amp; \sum_{v:(u,v)\in A}p_{uv}\leq \frac{1}{\ell}Q_u &amp; \forall u\in U\\ &amp; \frac{1}{\ell}\sum_{k\in K}q_kx_{uv}^k\leq p_{uv} &amp; \forall (u,v)\in A \end{align}</span></p> <p>Although these constraint work in theory, in practice they hinder the scalability of my model. Now I could simply drop the rounding operator and approximate the capacity constraints: <span class="math-container">\begin{align} &amp; \sum_{v:(u,v)\in A} \sum_{k\in K}q_kx_{uv}^k\leq Q_u &amp; \forall u\in U \end{align}</span> but this creates significant capacity deviations. Here's a simple numerical example. Imagine a node <span class="math-container">$u\in U$</span> with 10 arcs emanating from this node. When we evaluate the term <span class="math-container">$\sum_{k\in K}q_kx_{uv}^k$</span> for each of these 10 arcs, we find the values: <span class="math-container">$4.04,0.2,0.2,\dots,0.2$</span>. If we were to evaluate the left hand side of the capacity constraint with <span class="math-container">$\ell=0.5$</span>, we find: <span class="math-container">$4.5+0.5+0.5+\dots+0.5=9$</span>. Without the rounding operator, we would get <span class="math-container">$4.04+0.2+0.2+\dots+0.2=5.84$</span> which is a significant underestimation. This deviation becomes worse when the number of arcs and commodities increases or when <span class="math-container">$\ell$</span> is set to 1.</p> <p>Is there a better way to model these capacity constraints? This is an industrial application: I wouldn't mind to over or underestimate the exact capacity by some margin if this would improve scalability of the model.</p> https://or.stackexchange.com/q/7333 -3 No value for uninitialized NumericValue object Cindy Paola Guzman Lascano https://or.stackexchange.com/users/7581 2021-11-23T21:23:14Z 2021-11-26T21:33:31Z <p>I'm working on an optimization model in python with the pyomo library. However I'm getting an error message in python that I cannot seem to understand. The code and error message is below. My code is</p> <pre><code>from pyomo.environ import * from pyomo.opt import SolverFactory import json model = ConcreteModel() with open('C:/Users/cindy/python/Pyomo/ADMM_Distributed_Pyomo/input_data.json') as f: par = json.load(f) # Declaro os parametros model.P_PV = Param(initialize = par['P_PV']) model.P_D = Param(initialize = par['P_D']) model.P_SE = Param(initialize = par['P_SE']) model.P_ESS_max = Param(initialize = par['P_ESS_max']) model.E_ESS_max = Param(initialize = par['E_ESS_max']) model.P_load_max = Param(initialize = par['P_load_max']) model.E0 = Param(initialize = par['E0']) model.delta = Param(initialize = par['delta']) model.custo_venda = Param(initialize = par['custo_venda']) model.custo_compra = Param(initialize = par['custo_compra']) model.custo_load = Param(initialize = par['custo_load']) model.rho = Param(initialize = par['rho']) model.epsilon = Param(initialize = par['epsilon']) model.iter = Param(initialize = par['iter']) model.tol_lambda = Param(initialize = par['tol_lambda']) model.tol_var = Param(initialize = par['tol_var']) model.P_SE_in_1_param = Param(initialize = par['P_SE_in_1_param']) model.P_SE_out_1_param = Param(initialize = par['P_SE_out_1_param']) model.P_ESS_1_param = Param(initialize = par['P_ESS_1_param']) model.E_ESS_1_param = Param(initialize = par['E_ESS_1_param']) model.P_load_1_param = Param(initialize = par['P_load_1_param']) model.P_SE_in_2_param = Param(initialize = par['P_SE_in_2_param']) model.P_SE_out_2_param = Param(initialize = par['P_SE_out_2_param']) model.P_ESS_2_param = Param(initialize = par['P_ESS_2_param']) model.E_ESS_2_param = Param(initialize = par['E_ESS_2_param']) model.P_load_2_param = Param(initialize = par['P_load_2_param']) model.P_SE_in_3_param = Param(initialize = par['P_SE_in_3_param']) model.P_SE_out_3_param = Param(initialize = par['P_SE_out_3_param']) model.P_ESS_3_param = Param(initialize = par['P_ESS_3_param']) model.E_ESS_3_param = Param(initialize = par['E_ESS_3_param']) model.P_load_3_param = Param(initialize = par['P_load_3_param']) model.lambda_P_SE_in_a = Param(initialize = par['lambda_P_SE_in_a']) model.lambda_P_SE_out_a = Param(initialize = par['lambda_P_SE_out_a']) model.lambda_P_ESS_a = Param(initialize = par['lambda_P_ESS_a']) model.lambda_E_ESS_a = Param(initialize = par['lambda_E_ESS_a']) model.lambda_P_load_a = Param(initialize = par['lambda_P_load_a']) model.lambda_P_SE_in_b = Param(initialize = par['lambda_P_SE_in_b']) model.lambda_P_SE_out_b = Param(initialize = par['lambda_P_SE_out_b']) model.lambda_P_ESS_b = Param(initialize = par['lambda_P_ESS_b']) model.lambda_E_ESS_b = Param(initialize = par['lambda_E_ESS_b']) model.lambda_P_load_b = Param(initialize = par['lambda_P_load_b']) model.lambda_P_SE_in_a_ant = Param(initialize = par['lambda_P_SE_in_a_ant']) model.lambda_P_SE_out_a_ant = Param(initialize = par['lambda_P_SE_out_a_ant']) model.lambda_P_ESS_a_ant = Param(initialize = par['lambda_P_ESS_a_ant']) model.lambda_E_ESS_a_ant = Param(initialize = par['lambda_E_ESS_a_ant']) model.lambda_P_load_a_ant = Param(initialize = par['lambda_P_load_a_ant']) model.lambda_P_SE_in_b_ant = Param(initialize = par['lambda_P_SE_in_b_ant']) model.lambda_P_SE_out_b_ant = Param(initialize = par['lambda_P_SE_out_b_ant']) model.lambda_P_ESS_b_ant = Param(initialize = par['lambda_P_ESS_b_ant']) model.lambda_E_ESS_b_ant = Param(initialize = par['lambda_E_ESS_b_ant']) model.lambda_P_load_b_ant = Param(initialize = par['lambda_P_load_b_ant']) #Declaracao das variaveis model.P_SE_in_1 = Var(bounds = (0,None),within=NonNegativeReals) #Potência importada 1 model.P_SE_out_1 = Var(bounds = (0,None),within=NonNegativeReals) #Potência vendida 1 model.P_ESS_1 = Var(bounds = (None,None)) #Potência do BESS 1 model.E_ESS_1 = Var(bounds = (0,None),within=NonNegativeReals) #Energia do BESS 1 model.P_load_1 = Var(bounds = (0,None),within=NonNegativeReals) #Potência da carga controlada 1 model.P_SE_in_2 = Var(bounds = (0,None),within=NonNegativeReals) #Potência importada 2 model.P_SE_out_2 = Var(bounds = (0,None),within=NonNegativeReals)#Potência vendida 2 model.P_ESS_2 = Var(bounds = (None,None))#Potência do BESS 2 model.E_ESS_2 = Var(bounds = (0,None),within=NonNegativeReals)#Energia do BESS 2 model.P_load_2 = Var(bounds = (0,None),within=NonNegativeReals)#Potência da carga controlada 2 model.P_SE_in_3 = Var(bounds = (0,None),within=NonNegativeReals)#Potência importada 3 model.P_SE_out_3 = Var(bounds = (0,None),within=NonNegativeReals)#Potência vendida 3 model.P_ESS_3 = Var(bounds = (None,None))#Potência do BESS 3 model.E_ESS_3 = Var(bounds = (0,None),within=NonNegativeReals)#Energia do BESS 3 model.P_load_3 = Var(bounds = (0,None),within=NonNegativeReals)#Potência da carga controlada 3 #Declaro as funçoes objetivo #Declaro a FO1 model.obj1 = Objective(expr = model.delta * model.custo_compra * model.P_SE_in_1 - model.delta * model.custo_venda * model.P_SE_out_1 + model.delta * model.custo_load*(model.P_load_max - model.P_load_3_param) + model.lambda_P_SE_in_a*(model.P_SE_in_1 - model.P_SE_in_2_param) + model.rho/2*(model.P_SE_in_1 - model.P_SE_in_2_param)**2 + model.lambda_P_SE_out_a*(model.P_SE_out_1 - model.P_SE_out_2_param) + model.rho/2*(model.P_SE_out_1 - model.P_SE_out_2_param)**2 + model.lambda_P_ESS_a*(model.P_ESS_1 - model.P_ESS_2_param) + model.rho/2*(model.P_ESS_1 - model.P_ESS_2_param)**2 + model.lambda_E_ESS_a*(model.E_ESS_1 - model.E_ESS_2_param) + model.rho/2*(model.E_ESS_1 - model.E_ESS_2_param)**2 + model.lambda_P_load_a*(model.P_load_1 - model.P_load_2_param) + model.rho/2*(model.P_load_1 - model.P_load_2_param)**2 + model.lambda_P_SE_in_b*(model.P_SE_in_1 - model.P_SE_in_3_param) + model.rho/2*(model.P_SE_in_1 - model.P_SE_in_3_param)**2 + model.lambda_P_SE_out_b*(model.P_SE_out_1 - model.P_SE_out_3_param) + model.rho/2*(model.P_SE_out_1 - model.P_SE_out_3_param)**2 + model.lambda_P_ESS_b*(model.P_ESS_1 - model.P_ESS_3_param) + model.rho/2*(model.P_ESS_1 - model.P_ESS_3_param)**2 + model.lambda_E_ESS_b*(model.E_ESS_1 - model.E_ESS_3_param) + model.rho/2*(model.E_ESS_1 - model.E_ESS_3_param)**2 + model.lambda_P_load_b*(model.P_load_1 - model.P_load_3_param) + model.rho/2*(model.P_load_1 - model.P_load_3_param)**2) #Declaro as restricões da FO1 model.r1 = Constraint(expr = model.P_SE_in_1 + model.P_PV == model.P_ESS_1 + model.P_D + model.P_SE_out_1 + model.P_load_1) model.r2 = Constraint(expr = model.P_SE_in_1 &lt;= model.P_SE) model.r3 = Constraint(expr = model.P_SE_out_1 &lt;= model.P_SE) #Declaro a FO2 model.obj2 = Objective(expr = model.lambda_P_SE_in_a*(model.P_SE_in_1_param - model.P_SE_in_2) + model.rho/2*(model.P_SE_in_1_param - model.P_SE_in_2)**2 + model.lambda_P_SE_out_a*(model.P_SE_out_1_param - model.P_SE_out_2) + model.rho/2*(model.P_SE_out_1_param - model.P_SE_out_2)**2 + model.lambda_P_ESS_a*(model.P_ESS_1_param - model.P_ESS_2) + model.rho/2*(model.P_ESS_1_param - model.P_ESS_2)**2 + model.lambda_E_ESS_a*(model.E_ESS_1_param - model.E_ESS_2) + model.rho/2*(model.E_ESS_1_param - model.E_ESS_2)**2 + model.lambda_P_load_a*(model.P_load_1_param - model.P_load_2) + model.rho/2*(model.P_load_1_param - model.P_load_2)**2) #Declaro as restricões da FO2 model.r5 = Constraint(expr = model.E_ESS_2 == model.E0 + model.delta * model.P_ESS_2) model.r6 = Constraint(expr = -1*model.P_ESS_max &lt;= model.P_ESS_2) model.r7 = Constraint(expr = model.P_ESS_2 &lt;= model.P_ESS_max) model.r8 = Constraint(expr = model.E_ESS_2 &lt;= model.E_ESS_max) #Declaro a FO3 model.obj3 = Objective(expr = model.delta * model.custo_load *(model.P_load_max - model.P_load_3) + model.lambda_P_SE_in_b*(model.P_SE_in_1_param - model.P_SE_in_3) + model.rho/2*(model.P_SE_in_1_param - model.P_SE_in_3)**2 + model.lambda_P_SE_out_b*(model.P_SE_out_1_param - model.P_SE_out_3) +model.rho/2*(model.P_SE_out_1_param - model.P_SE_out_3)**2 + model.lambda_P_ESS_b*(model.P_ESS_1_param - model.P_ESS_3) + model.rho/2*(model.P_ESS_1_param - model.P_ESS_3)**2 + model.lambda_E_ESS_b*(model.E_ESS_1_param - model.E_ESS_3) + model.rho/2*(model.E_ESS_1_param - model.E_ESS_3)**2 + model.lambda_P_load_b*(model.P_load_1_param - model.P_load_3) + model.rho/2*(model.P_load_1_param - model.P_load_3)**2) #Declaro as restricões da FO3 model.r4 = Constraint(expr = model.P_load_3 &lt;= model.P_load_max) #Problem solution while model.tol_lambda &gt; 10e5 and model.tol_var &gt; 10e5 : model.lambda_P_SE_in_a_ant = model.lambda_P_SE_in_a model.lambda_P_SE_out_a_ant = model.lambda_P_SE_out_a model.lambda_P_ESS_a_ant = model.lambda_P_ESS_a model.lambda_E_ESS_a_ant = model.lambda_E_ESS_a model.lambda_P_load_a_ant = model.lambda_P_load_a model.lambda_P_SE_in_b_ant = model.lambda_P_SE_in_b model.lambda_P_SE_out_b_ant = model.lambda_P_SE_out_b model.lambda_P_ESS_b_ant = model.lambda_P_ESS_b model.lambda_E_ESS_b_ant = model.lambda_E_ESS_b model.lambda_P_load_b_ant = model.lambda_P_load_b #Solucao da FO1 opt = SolverFactory('glpk') model.obj1.activate() model.obj2.deactivate() model.obj3.deactivate() results = opt.solve(model) model.P_SE_in_1_param == model.P_SE_in_1 model.P_SE_out_1_param == model.P_SE_out_1 model.P_ESS_1_param == model.P_ESS_1 model.E_ESS_1_param == model.E_ESS_1 model.P_load_1_param == model.P_load_1 #Solucao da FO2 ''' opt = SolverFactory('glpk') model.obj1.deactivate() model.obj3.deactivate() results = opt.solve(model) ''' model.P_SE_in_2_param == model.P_SE_in_2 model.P_SE_out_2_param == model.P_SE_out_2 model.P_ESS_2_param == model.P_ESS_2 model.E_ESS_2_param == model.E_ESS_2 model.P_load_2_param == model.P_load_2 #Solucao da FO3 ''' opt = SolverFactory('glpk') model.obj1.deactivate() model.obj2.deactivate() results = opt.solve(model) ''' model.P_SE_in_3_param == model.P_SE_in_3 model.P_SE_out_3_param == model.P_SE_out_3 model.P_ESS_3_param == model.P_ESS_3 model.E_ESS_3_param == model.E_ESS_3 model.P_load_3_param == model.P_load_3 #atualizacao variavel dual model.lambda_P_SE_in_a = model.lambda_P_SE_in_a_ant + model.rho * (model.P_SE_in_1_param - model.P_SE_in_2_param) model.lambda_P_SE_out_a = model.lambda_P_SE_out_a_ant + model.rho * (model.P_SE_out_1_param - model.P_SE_out_2_param) model.lambda_P_ESS_a = model.lambda_P_ESS_a_ant + model.rho * (model.P_ESS_1_param - model.P_ESS_2_param) model.lambda_E_ESS_a = model.lambda_E_ESS_a_ant + model.rho * (model.E_ESS_1_param - model.E_ESS_2_param) model.lambda_P_load_a = model.lambda_P_load_a_ant + model.rho * (model.P_load_1_param - model.P_load_2_param) model.lambda_P_SE_in_b = model.lambda_P_SE_in_b_ant + model.rho * (model.P_SE_in_1_param - model.P_SE_in_3_param) model.lambda_P_SE_out_b = model.lambda_P_SE_out_b_ant + model.rho * (model.P_SE_out_1_param - model.P_SE_out_3_param) model.lambda_P_ESS_b = model.lambda_P_ESS_b_ant + model.rho * (model.P_ESS_1_param - model.P_ESS_3_param) model.lambda_E_ESS_b = model.lambda_E_ESS_b_ant + model.rho * (model.E_ESS_1_param - model.E_ESS_3_param) model.lambda_P_load_b = model.lambda_P_load_b_ant + model.rho * (model.P_load_1_param - model.P_load_3_param) # Calcula os critérios de parada model.tol_var = abs(model.P_SE_in_1_param - model.P_SE_in_2_param) + abs(model.P_SE_out_1_param - model.P_SE_out_2_param) + abs(model.P_ESS_1_param - model.P_ESS_2_param) + abs(model.E_ESS_1_param - model.E_ESS_2_param) + abs(model.P_load_1_param - model.P_load_2_param) + abs(model.P_SE_in_1_param - model.P_SE_in_3_param) + abs(model.P_SE_out_1_param - model.P_SE_out_3_param) + abs(model.P_ESS_1_param - model.P_ESS_3_param) + abs(model.E_ESS_1_param - model.E_ESS_3_param) + abs(model.P_load_1_param - model.P_load_3_param) model.tol_lambda = abs(model.lambda_P_SE_in_a - model.lambda_P_SE_in_a_ant) + abs(model.lambda_P_SE_out_a - model.lambda_P_SE_out_a_ant) + abs(model.lambda_P_ESS_a - model.lambda_P_ESS_a_ant) + abs(model.lambda_E_ESS_a - model.lambda_E_ESS_a_ant) + abs(model.lambda_P_load_a - model.lambda_P_load_a_ant) + abs(model.lambda_P_SE_in_b - model.lambda_P_SE_in_b_ant) + abs(model.lambda_P_SE_out_b - model.lambda_P_SE_out_b_ant) + abs(model.lambda_P_ESS_b - model.lambda_P_ESS_b_ant) + abs(model.lambda_E_ESS_b - model.lambda_E_ESS_b_ant) + abs(model.lambda_P_load_b - model.lambda_P_load_b_ant) iter = iter + 1, </code></pre> <p>The error is</p> <pre><code>ERROR: evaluating object as numeric value: P_SE_in_1 (object: &lt;class 'pyomo.core.base.var.ScalarVar'&gt;) No value for uninitialized NumericValue object P_SE_in_1 Traceback (most recent call last): File &quot;C:\Users\cindy\python\Pyomo\ADMM_Distributed_Pyomo\ADMM_Distributed_Pyomo.py&quot;, line 218, in &lt;module&gt; FO1 = model.obj1.expr() File &quot;pyomo\core\expr\numeric_expr.pyx&quot;, line 218, in pyomo.core.expr.numeric_expr.ExpressionBase.__call__ File &quot;C:\Users\cindy\AppData\Local\Packages\PythonSoftwareFoundation.Python.3.9_qbz5n2kfra8p0\LocalCache\local-packages\Python39\site-packages\pyomo\core\expr\visitor.py&quot;, line 1045, in evaluate_expression return visitor.dfs_postorder_stack(exp) File &quot;C:\Users\cindy\AppData\Local\Packages\PythonSoftwareFoundation.Python.3.9_qbz5n2kfra8p0\LocalCache\local-packages\Python39\site-packages\pyomo\core\expr\visitor.py&quot;, line 572, in dfs_postorder_stack flag, value = self.visiting_potential_leaf(_sub) File &quot;C:\Users\cindy\AppData\Local\Packages\PythonSoftwareFoundation.Python.3.9_qbz5n2kfra8p0\LocalCache\local-packages\Python39\site-packages\pyomo\core\expr\visitor.py&quot;, line 953, in visiting_potential_leaf return True, value(node, exception=self.exception) File &quot;pyomo\core\expr\numvalue.pyx&quot;, line 156, in pyomo.core.expr.numvalue.value File &quot;pyomo\core\expr\numvalue.pyx&quot;, line 143, in pyomo.core.expr.numvalue.value ValueError: No value for uninitialized NumericValue object P_SE_in_1 PS C:\Users\cindy\python\Pyomo\ADMM_Distributed_Pyomo&gt; </code></pre> https://or.stackexchange.com/q/7332 1 Importance of Polynomial Time Approximation algorithm for special case Shawky https://or.stackexchange.com/users/7538 2021-11-23T20:27:04Z 2021-11-23T20:27:04Z <p>I am an MSc student, searching for an interesting topic in operations research for NP-hard problems. I need your opinion and recommendation about: &quot;Searching a polynomial time approximation algorithm for the special case of facility layout problem: single row equidistant facility layout.&quot;</p> <p>Regarding the literature, there is no polynomial time approximation algorithm for this problem, only heuristic, metaheuristic, and exact algorithms (solve up to n=35). My inquires are:</p> <p>Is this topic interesting in scientific research? What are the most expectations of results for such an algorithm to be acceptable for publication? Any recommendations for journals names in this scope?</p> https://or.stackexchange.com/q/7331 4 Chance constrained optimization - interpretation Djames https://or.stackexchange.com/users/609 2021-11-22T14:02:54Z 2021-11-23T05:26:05Z <p>Suppose that we have a stochastic vector <span class="math-container">$\psi$</span> and <span class="math-container">$S$</span> realisations of <span class="math-container">$\psi$</span> given by <span class="math-container">$\psi_1,\dots,\psi_S$</span> with equal probability of occurrence. In addition, we have constraints of the form <span class="math-container">\begin{equation} h_i(x,\psi)\leq b_i,\quad \forall i=1,...,m \end{equation}</span> for the decision vector <span class="math-container">$x$</span>.</p> <p>A <em>joint chance constraint</em> is then given by <span class="math-container">\begin{equation} P(\ h_i(x,\psi)\leq b_i , \quad \forall i=1,..,m\ )\geq \alpha \end{equation}</span> stating that we can accept that some (or all) of these constraints are violation with a probability of <span class="math-container">$1-\alpha$</span>. We could also write <em>single chance constraints</em> as follows <span class="math-container">\begin{equation} P(\ h_i(x,\psi)\leq b_i \ )\geq \alpha, \quad \forall i=1,..,m \end{equation}</span> stating that we will accept violations of the individual constraints with a probability of <span class="math-container">$1-\alpha$</span>. Using binary variables <span class="math-container">$z^s$</span> equalling 0 <em>iff</em> all constraints are satisfied in realisation <span class="math-container">$s$</span>, we can formulate the joint chance constraint as the MIP <span class="math-container">\begin{align} h_i(x,\psi_s)\leq b_i+Mz^s,&amp;&amp;\forall i=1,...,m,s=1,...,S\\ \sum_{s=1}^Sz^s\leq \lfloor (1-\alpha)S\rfloor \end{align}</span></p> <p>Using binary variables <span class="math-container">$z^s_i$</span> equalling 0 <em>iff</em> constraint <span class="math-container">$i$</span> is statisfied in realisation <span class="math-container">$s$</span> we can formulate the single chance constraint version as follows: <span class="math-container">\begin{align} &amp;h_i(x,\psi_s)\leq b_i+Mz^s_i,&amp;&amp;\forall i=1,...,m,s=1,...,S\\ &amp;\sum_{s=1}^Sz^s_i\leq \lfloor (1-\alpha)S\rfloor,&amp;&amp;\forall i=1,...,m \end{align}</span></p> <p>My question is, what is the interpretation of the following MIP <span class="math-container">\begin{align} &amp;h_i(x,\psi_s)\leq b_i+Mz^s_i,&amp;&amp;\forall i=1,...,m,s=1,...,S\\ &amp;\sum_{i=1}^m\sum_{s=1}^Sz^s_i\leq \lfloor (1-\alpha)S\rfloor \end{align}</span> Does it have some sensible interpretation?</p> https://or.stackexchange.com/q/7330 1 Sources of Min-Cost Flow Models That Utilize Binary Variables for Transportation Networks Mert Saner https://or.stackexchange.com/users/7526 2021-11-22T12:06:26Z 2021-11-22T12:13:19Z <p>I am looking for articles that include min-cost flow models with binary variables for flow transportation like gas networks, traffic systems, heating systems. Is there any specific place(like OR Library) where I can find all related articles? Can you share articles you know, that use min-cost flow models which implement binary variables for transportation networks?</p> https://or.stackexchange.com/q/7327 4 Python "Coffee Shop Scheduling Problem" - Scheduling Lunches/Breaks t25 https://or.stackexchange.com/users/7573 2021-11-21T20:32:44Z 2021-11-21T22:24:19Z <p>I'm working on an employee scheduling program in python. Having never done this before, I've been researching different libraries that can be used to accomplish the task.</p> <p>Unfortunately, none of the examples I have been able to find for scheduling problems (for any library) seem to address one of the most common scheduling requirements - optimally scheduling breaks and lunches.</p> <p>My program is ultimately going to fit the &quot;coffee shop scheduling&quot; model found here: <a href="https://towardsdatascience.com/how-to-solve-a-staff-scheduling-problem-with-python-63ae50435ba4" rel="nofollow noreferrer">https://towardsdatascience.com/how-to-solve-a-staff-scheduling-problem-with-python-63ae50435ba4</a></p> <p>I know how many employees I need to meet demand every hour, but it swings wildly. So I need the program to select the optimal shifts to staff, from a list of known shifts, to be able to meet all demand.</p> <p>For reference, number of employees needed would look something like this:</p> <p><a href="https://i.stack.imgur.com/AxdoN.png" rel="nofollow noreferrer"><img src="https://i.stack.imgur.com/AxdoN.png" alt="enter image description here" /></a></p> <p>etc.</p> <p>I'd have a list of shifts to pull from, with pre-defined breaks. It's essentially a cartesian product of all acceptable shift/break/lunch combinations:</p> <p><a href="https://i.stack.imgur.com/X1bpG.png" rel="nofollow noreferrer"><img src="https://i.stack.imgur.com/X1bpG.png" alt="enter image description here" /></a></p> <p>etc.,</p> <p>I'm sure there are a number of ways to do this, but what would be one possible way to incorporate the break/lunch schedules into PuLP, Pyomo, or Google OR-Tools? I've researched PuLP, Pyomo, and Google OR-Tools and I'm willing to use whichever library can best solve for this problem, they all seem like they can do the rest of what I need.</p> <p>Any help or general direction would be greatly appreciated. I'm a little bit stuck at this juncture.</p> https://or.stackexchange.com/q/7324 6 How to keep solutions stable/persistent in a problem with many equally good solutions? Geoffrey Brent https://or.stackexchange.com/users/546 2021-11-21T02:42:39Z 2021-11-23T23:37:22Z <p>Suppose we have a worker-assignment problem where we seek to assign Alice, Bob, Chris, ... to jobs A, B, C, ... subject to various constraints and some objective function based on these assignments, and with the following considerations:</p> <ol> <li>Many employees are identical except in name. (Workers <span class="math-container">$i$</span> and <span class="math-container">$j$</span> considered &quot;identical&quot; if for every feasible solution, the allocations for worker <span class="math-container">$i$</span> can be swapped with those for <span class="math-container">$j$</span> without changing the feasibility or objective-function value of the solution).</li> <li>Similarly, many jobs are identical except in name.</li> </ol> <p>Hence, there will be many optimal solutions that differ from one another only by swaps between identical entities - let's call these &quot;sibling&quot; solutions.</p> <ol start="3"> <li><p><strike>Stability/reproducibility</strike> Persistence of solutions is desired. If we rerun the same problem to optimality, we want to get the exact same solution, not one of its siblings. If we make small changes to the problem and rerun, we want to change the solution by enough to achieve optimality in the new version of the problem, but no more than that.</p> </li> <li><p>We may wish to do this many times, without keeping records of all the old solutions.</p> </li> </ol> <p>How can we achieve this?</p> <p>If it weren't for #4, one possible solution would be to add a small term to the objective function that penalises changes from previous solutions. Unfortunately this becomes impractical if we are doing a lot of runs with small changes.</p> <p>I have figured out one way to do this and will post as a self-answer, but I am interested in hearing others.</p> https://or.stackexchange.com/q/7321 2 How to normalize the objective functions of multi-objective optimization into uniform form? MAJID majid https://or.stackexchange.com/users/7281 2021-11-20T11:21:18Z 2021-11-20T17:14:01Z <p>In my bi-objective model, the range of solution value for the first objective is large than the second objective. I decide to obtain a single solution by the weighted sum approach and solve it using metaheuristics. How can I do that?</p> https://or.stackexchange.com/q/7319 3 Inverse of weighted sum of positive definite matrices Jagadish https://or.stackexchange.com/users/7561 2021-11-19T20:55:31Z 2021-11-20T09:15:45Z <p>Let us suppose <span class="math-container">$I_1, \ldots, I_n$</span> are symmetric and positive definite matrices. Let <span class="math-container">$\mathbf{u}$</span> be the vector with <span class="math-container">$n$</span> 1s. I'm interested in the following optimization problem:</p> <p><span class="math-container">$$\min \; u^T (x_1I_1+x_2I_2+\ldots + x_nI_n)^{-1}u$$</span></p> <p>such that <span class="math-container">$0 \leq x_i\leq 1$</span> and <span class="math-container">$\sum x_i \leq 1$</span>.</p> <p><em>Question:</em> Is the above optimization problem convex over the domain of <span class="math-container">$x_i$</span>s?</p> https://or.stackexchange.com/q/7318 3 Multiple SKU forecast for Intermittent Demand Dan https://or.stackexchange.com/users/7558 2021-11-19T11:37:13Z 2021-11-25T10:20:28Z <p>I've been tasked to generate a forecast for our newly operation business which has more than 500+ sku. Almost 90% of them are following intermittent demand pattern, with very few data points to train upon (maximum I recorded 80-90 data points for each product) I cannot train an ML model.</p> <p>Is there a way or model where I can train such small data set using excel or any other tool?</p> <p>Link to data <a href="https://docs.google.com/spreadsheets/d/1W-wmYGzLtZpqV0IcYpeV6MrgdKESGmHr/edit?usp=sharing&amp;ouid=104686423957654811582&amp;rtpof=true&amp;sd=true" rel="nofollow noreferrer">Intermittent_DATA</a></p> https://or.stackexchange.com/q/7315 5 Optimize for bonuses within a group (knapsack) Eddie https://or.stackexchange.com/users/7539 2021-11-19T04:09:14Z 2021-11-25T16:54:26Z <p>I am trying to create an LP problem which is like the knapsack problem but with groups of items. Let's say there are 10 groups of items each with up to 5 items. Each group has one special item and you must choose only one, let's call that group the &quot;special group&quot;. There's a bonus of 4 if you take 1 more item from the special group and a bonus of 8 if you pick two.</p> <p>What I have tried is having <code>g[1:10]</code> which identifies the special group, then <code>a[5row, 10col]</code> which tells me how many items were taken from each group, then a new variable <code>b[1:5]</code> (max # of items). I then tried this constraint: <code>[i=1:5], b[i] &lt;= sum(g[j] * a[i,j] for j in 1:10)</code>. My goal here was that since I only care about my special group, I would multiply each group's item count in <code>a</code> by whether it's the special group in <code>g</code>. This introduces a quadratic constraint which is not allowed in my solver.</p> <p>What is the right approach to this problem?</p> https://or.stackexchange.com/q/7313 2 How to model a cumulative resource constraint with same family condition? LouisPopovic https://or.stackexchange.com/users/7347 2021-11-18T19:18:53Z 2021-11-19T13:18:39Z <p>Assume that we have 3 tasks to schedule : <span class="math-container">$x_{a1}, x_{a2}, x_{b3}$</span>. They all use the same cumulative resource : <span class="math-container">$r_1$</span>. Each activity increases <span class="math-container">$r_1$</span> by <span class="math-container">$1$</span> at it starts time and decreases <span class="math-container">$r_1$</span> by <span class="math-container">$1$</span> at it's end time. <span class="math-container">$r_1$</span> has a maximum capacity of <span class="math-container">$1$</span>. Where things get complicated is that <span class="math-container">$x_{a1}$</span> and <span class="math-container">$x_{a2}$</span> can be executed at the same time and only consume <span class="math-container">$1$</span> resource of <span class="math-container">$r_1$</span>. So in this example, <span class="math-container">$x_{a1}$</span> and <span class="math-container">$x_{a2}$</span> can be executed at the same time while <span class="math-container">$x_{b3}$</span> must be executed before of after both <span class="math-container">$x_{a1}$</span> and <span class="math-container">$x_{a2}$</span> are executed. Basically, tasks from the same family only consume <span class="math-container">$1$</span> resource even if multiple tasks from this family are overlapping.</p> <p>Note that the maximum capacity of the ressource could be <span class="math-container">$&gt;1$</span></p> <p>With CPOptimizer I know how to use <code>cumulFunction()</code> but I don't know how I could implement the fact that multiple activites can be scheduled simultaneously while only consuming one resource.</p> https://or.stackexchange.com/q/7308 4 Is it possible to identify all possible Irreducible Infeasible Sets (IIS) for an infeasible Integer Linear Programming problem? (ILP)? Ramy Fouad https://or.stackexchange.com/users/7554 2021-11-18T03:37:28Z 2021-11-18T13:58:35Z <p>For an Integer Linear Programming problem (ILP), an irreducible infeasible set (IIS) is an infeasible subset of constraints, variable bounds, and integer restrictions that becomes feasible if any single constraint, variable bound, or integer restriction is removed. It is possible to have more than one IIS in an infeasible ILP. Is it possible to identify all possible Irreducible Infeasible Sets (IIS) for an infeasible Integer Linear Programming problem (ILP)? Ideally, I aim to find the MIN IIS COVER, which is the smallest cardinality subset of constraints to remove such that at least one constraint is removed from every IIS.</p> https://or.stackexchange.com/q/7303 1 Capacitated Facility Location problem with a dynamic set up cost Shibaprasadb https://or.stackexchange.com/users/5532 2021-11-17T14:54:19Z 2021-11-18T16:22:58Z <p>I am trying to solve a Capacitated Facility Location Problem (CFLP) with a dynamic setup cost in <code>R</code>.</p> <p>The problem statement is this:</p> <ol> <li>I have the transport cost</li> <li>The fixed operating cost (manual labor and other expenses) is known</li> <li>I know the dropping points with loads and all the details</li> <li>The per square ft. cost of rent of a place is known</li> <li>The size of the Facility will be a function of the load. So the rent will depend on how much load is getting allocated in that place.</li> </ol> <p>Assuming the rent will vary like this:</p> <p>rent= rent_per_square_ft * load* 0.10</p> <p>Now, I have accommodated the first 4 conditions in my code. But I am not sure how the number 5 can be accommodated.</p> <p>My model looks like this in <code>R</code>(if it can be of any help):</p> <pre><code>#m is the number of potential facility/service center (SC) locations #n is the number of customer locations model &lt;- MIPModel() %&gt;% # 1 if customer i gets assigned to SC j add_variable(x[i, j], i = 1:n, j = 1:m, type = &quot;binary&quot;) %&gt;% # 1 if SC j is built add_variable(y[j], j = 1:m, type = &quot;binary&quot;) %&gt;% # Objective function set_objective(sum_expr(transportcost(i, j) * x[i, j], i = 1:n, j = 1:m) + sum_expr(fixedcost[j] * y[j], j = 1:m), &quot;min&quot;) %&gt;% # Every customer needs to be assigned to a SC add_constraint(sum_expr(x[i, j], j = 1:m) == 1, i = 1:n) %&gt;% # If a customer is assigned to a SC, then the SC must be built add_constraint(x[i,j] &lt;= y[j], i = 1:n, j = 1:m) %&gt;% #The demand of customers shouldn't exceed SC capacities add_constraint(sum_expr(demand[i] * x[i, j], i = 1:n) &lt;= capacity[j] * y[j], j = 1:m) </code></pre> <p>I am looking for any headway. Even any link to a relevant article might help.</p> https://or.stackexchange.com/q/7302 3 Conditional constraint with a strict inequality athing https://or.stackexchange.com/users/7544 2021-11-17T10:30:03Z 2021-11-18T13:31:21Z <p>It's almost this question: <a href="https://or.stackexchange.com/questions/4715/formulating-the-conditional-constraint">Formulating the conditional constraint</a></p> <p>But there they have non-strict inequality. I have <span class="math-container">$x_i$</span> a boolean decision var and <span class="math-container">$Q_i$</span> as a nonnegative integer decision variable such that</p> <ol> <li>if <span class="math-container">$x_i = 0$</span>, then <span class="math-container">$Q_i = 0$</span></li> <li>if <span class="math-container">$x_i = 1$</span>, then <span class="math-container">$Q_i \gt 0$</span> (note the strict inequality!).</li> </ol> <p>Lets say I dont have upper bound on <span class="math-container">$Q_i$</span>; is there a mathematical relation between <span class="math-container">$x_i$</span> and <span class="math-container">$Q_i$</span> you can write directly or is the following way to go?</p> <pre><code>dvar boolean x[I]; dvar int+ Q[I]; subject to { forall(i in I) { (x[i]==0) =&gt; (Q[i] == 0); (x[i]==1) =&gt; (Q[i] &gt; 0); } </code></pre> <p>I can formulate the constraint other way but does it help I'm not sure:</p> <ol> <li>if <span class="math-container">$Q_i = 0$</span>, then <span class="math-container">$x_i = 0$</span></li> <li>if <span class="math-container">$Q_i \gt 0$</span>, then <span class="math-container">$x_i = 1$</span>.</li> </ol> https://or.stackexchange.com/q/7301 3 How to formulate this optimisation problem that has time varying components? asett https://or.stackexchange.com/users/7543 2021-11-17T08:36:58Z 2021-11-20T14:53:54Z <p>I've been working on a optimisation problem at work regarding blending. Essentially we have different types of grain stock with different properties and we bought them for different amounts. We then have a bunch of customer contracts requesting some amount of grain with certain property constraints (e.g. protein content must be less than this amount). We can blend our stocks to fulfill a contract provided that the property constraints of the contract are not violated. We can formulate this pretty easily as an LP problem where our objective function is to minimise the total price of grain stock we allocated to all of our contracts, and the constraints are given by the contracts property constraints, the grain amount specified by the contract and the amount of grain stock we have on hand.</p> <p>In reality however we are constantly receiving new stock throughout a given year and shipments for contracts will happen at different times. We will also receive new contracts throughout the year too. So at any point in time we potentially don't know stock and contracts we could have in the future. How would I take into account this time varying aspect to formulate this as an optimisation problem such that we could maximise profit for a given year? We were considering just optimising for each batch of contracts that correspond to the same shipping period individually using stock we would have at that point in time but is there a way to formulate this problem as one single global problem?</p> https://or.stackexchange.com/q/7285 1 Rule of thumb for introducing known constraints in Mixed Integer problem when they grow the formula worldsmithhelper https://or.stackexchange.com/users/4777 2021-11-15T21:09:31Z 2021-11-18T10:09:35Z <p>When encoding a problem into a Mixed Integer formula one faces a trade-off between embedding domain knowledge which might require new helper variables to express the constraints thereby growing the formulation.</p> <p><strong>Are you aware of any rule of thumb for mixed integer (Linear, Quadratic or Non-linear) problems whether encoding some additional constraint is worth the number of variables and constraints?</strong></p> <p>For MILP i might imagine it to look like the product of the simplex iteration time of the two variants times some function of the size of the possibilty space on both variants.</p> https://or.stackexchange.com/q/7113 7 How to model a non-overlap constraint between 2 groups of tasks? LouisPopovic https://or.stackexchange.com/users/7347 2021-10-14T19:58:41Z 2021-11-19T13:19:16Z <p>Let <span class="math-container">$T$</span> be a set of tasks. Each task <span class="math-container">$t \in T$</span> has a duration <span class="math-container">$d_t$</span>.</p> <p>Let <span class="math-container">$T_1 \subset T$</span> and <span class="math-container">$T_2 \subset T$</span> such as <span class="math-container">$T_1 \cap T_2 = \emptyset$</span></p> <p>How can I model the following constraint : all tasks from <span class="math-container">$T_1$</span> cannot overlap all tasks of <span class="math-container">$T_2$</span> and vice versa. So basically, tasks from a same set can overlap but when a task from a set is performed, no task from the other set can overlap it.</p> https://or.stackexchange.com/q/6602 1 How do I write the following constraints in Pyomo? user4387 https://or.stackexchange.com/users/0 2021-07-20T10:09:48Z 2021-11-19T13:21:15Z <p>I'm currently trying to replicate a vehicle routing problem paper, and I'm facing trouble writing the said constraints in the paper in Pyomo. I've been able to generate the data as required by the model, but unfortunately due to my limited experience with Pyomo, I'm unable to write some constraints.</p> <pre><code>I have declared my variables and index sets in Pyomo like this: J = list(S_i_dict.keys()) # Set of locations I = list(S_i_dict.keys()) # Set of Facilties K = [1,2,3,4,5,6,7,8,9,10,11,12,13,14,15] #Step 2: Define variables model_1.x_ij = Var(I,J, within = Binary) # if a demand point i is assigned to a vaccination centre j at a time period t model_1.y_i = Var(I, within = Binary) model_1.z_ijk = Var(I, J, K within = Binary) model_1.u_ik = Var(I, K, domain = NonNegativeIntegers) model_1.w_i = Var(I, domain = NonNegativeIntegers) </code></pre> <p><span class="math-container">$z_{iik} = 0 \ \ \forall i,k$</span></p> <p>I'm unable to write the above constraint as I've fixed <code>model_1.z_ijk = Var(I, J, K within = Binary)</code>.</p> <p><span class="math-container">$$U_{ik} - U_{jk} + pZ_{ijk} \leq p - W_{j} \forall \ \ i,j,k$$</span></p> <p>Similar to the above problem, I've fixed <code>model_1.u_ik = Var(I, K, domain = NonNegativeIntegers)</code>.</p> <p>Additionally, I'm finding trouble to write the following constraint -</p> <p><span class="math-container">$$W_i = \sum_{0 \leq j \leq n=1} b_j X_{ij} \forall \ \ i \leq n$$</span></p> <p>Hoping I receive some help on writing these constraints, I've already the Pyomo cookbook, but I continue to struggle.</p> https://or.stackexchange.com/q/6593 2 Consistency Level of Constraints Dav https://or.stackexchange.com/users/3418 2021-07-19T16:02:57Z 2021-11-19T13:20:57Z <p>Recently I started using CP and translated my first MIP model into a CP model with the IBM CP Optimizer. Lately I was confronted with a question I could not answer:</p> <p>&quot;What is the consistency level (e.g. bound consistent or domain consistent) applied to the constraints of the model?&quot;</p> <p>My Questions:</p> <ul> <li>What exactly is the consistency level of the constraints of the model?</li> <li>How can I find out about the consistency level of the constraints?</li> <li>For which reason would you want to let the reader know about the consistency level? or: Why would the reader be interested?</li> </ul> https://or.stackexchange.com/q/6402 4 Size timesteps for time-dependent TSP Theodeo https://or.stackexchange.com/users/5571 2021-06-07T08:25:08Z 2021-11-20T14:54:22Z <p>I'm solving a time-dependent travelling salesman problem. I have data from a whole year. I'm thinking about the size of the time steps. I think that the time-steps shouldn't be too large to be realistic (in 4 hours the traffic can change significantly) but I also think that it shouldn't be too small so that the data does not have the tendency to be scewed. Is there a &quot;best practice&quot; size for the TD-TSP or is smaller always better?</p> https://or.stackexchange.com/q/5884 2 Add Variables in Bulk for ortools Zachzhao https://or.stackexchange.com/users/5113 2021-03-09T23:51:12Z 2021-11-20T14:56:05Z <p>I am using the Python <code>ortools</code> to solve large LPs where the number of variables and constraints are in the hundred thousands (<span class="math-container">$10^5$</span>). Surprisingly, it takes longer to initialize variables and declare constraints than to solve the actual LP:</p> <div class="s-table-container"> <table class="s-table"> <thead> <tr> <th style="text-align: left;">Phase</th> <th>Average Time (secs)</th> </tr> </thead> <tbody> <tr> <td style="text-align: left;">Variable Initialization</td> <td>56.28</td> </tr> <tr> <td style="text-align: left;">Constraint Initialization</td> <td>203.73</td> </tr> <tr> <td style="text-align: left;">LP Solve</td> <td>116.34</td> </tr> </tbody> </table> </div> <p>Is there a way to lower the time to initialize variables and constraints by batching them? For my experiments, I solve many of these LPs, so any reduction in runtime would be appreciated.</p> https://or.stackexchange.com/q/5589 4 Safety Stock with Fill Rate Criterion TTY https://or.stackexchange.com/users/4743 2021-01-24T04:33:13Z 2021-11-22T10:39:55Z <p>When applying the <a href="https://en.wikipedia.org/wiki/Base_stock_model" rel="nofollow noreferrer">base stock inventory policy</a>, assuming the daily demands are normally distributed with parameter <span class="math-container">$(\mu, \sigma)$</span>, we can find the optimal parameter <span class="math-container">$S$</span> (the base stock level) in several different ways: (say we have 0 lead time and review period is 1 unit period, and assume infinite horizon)</p> <ol> <li>From the holding cost/stockout cost criterion: if holding cost per item is <span class="math-container">$h$</span> and stockout cost per item is <span class="math-container">$p$</span>, then <span class="math-container">$S = \mu + \sigma\Phi^{-1}(\frac{p}{p+h})$</span>, where <span class="math-container">$\Phi$</span> is the cdf of the standard normal distribution.</li> <li>If the stockout cost <span class="math-container">$p$</span> is difficult to estimate for the firm, then a service-level-based approach is used, in particular, the two most basic types of service levels are Cycle Service Level (type 1 service rate) and Fill Rate (type 2 service rate): To achieve a type 1 service level of <span class="math-container">$\alpha$</span>, we simply have <span class="math-container">$\alpha = \Phi((S-\mu)/\sigma)$</span>, so the base stock level <span class="math-container">$S = \mu + \sigma\Phi^{-1}(\alpha)$</span>.</li> <li>For a type 2 service level, the calculation is more complicated. The usual formula for approximating the fill rate is <span class="math-container">$\beta = 1-\frac{n(S)}{\mu}$</span>, where <span class="math-container">$n(S) = \sigma \mathcal{L}(z), z=(S-\mu)/\sigma$</span>, and <span class="math-container">$\mathcal{L}(z)$</span> is the standard normal loss function (<a href="https://books.google.com/books/about/Fundamentals_of_Supply_Chain_Theory.html?id=U7GTrLyVnPMC" rel="nofollow noreferrer">see e.g. in the appendix of this book</a>).</li> </ol> <p>The first two approaches give the base stock level a nice structure: <span class="math-container">$\mu$</span> is the cycle stock to cover the average demand in lead time and review period, while <span class="math-container">$\sigma \Phi^{-1}(\alpha)$</span> is the safety stock to buffer the fluctuations in lead time demand, and we have a rather simple description of the relation between the safety stock and the service level <span class="math-container">$\alpha$</span>. However, when using the fill rate approach, the base stock level <span class="math-container">$S$</span> is found by solving the nonlinear equation <span class="math-container">$\beta = 1-n(S)/\mu$</span>, we can still compute the associated base safety stock level <span class="math-container">$ss = S-\mu$</span>, but we no longer have a simple description on the relationship bewteen the safety stock level, the fill rate <span class="math-container">$\beta$</span> and the demand standard deviation <span class="math-container">$\sigma$</span>.</p> <p>So my question is: is there any approximate formula/asymptotic expression (as <span class="math-container">$\beta$</span> approaches <span class="math-container">$1$</span>) that gives the rough relation between the safety stock level <span class="math-container">$ss = S-\mu$</span> and the fill rate <span class="math-container">$\beta$</span> when <span class="math-container">$S$</span> is found by solving the fill rate constraint?</p> https://or.stackexchange.com/q/4485 6 How to simulate computational execution time? Matheus Diógenes Andrade https://or.stackexchange.com/users/2862 2020-07-05T17:32:48Z 2021-11-20T14:56:21Z <p>I am working with some computational experiments with an Integer Programming (IP) formulation over a well-known set of instances from the literature of my problem. And I would like to compare my formulation with another author formulation from the literature (making a benchmarking).</p> <p>However, the machine configuration in which the author's formulation was executed is different from mine. So, I thought, if I can find the equivalent time from the author's machine to mine, then I can use this equivalent time on my computational experiments, e.g if one hour of execution on the author's machine is equivalent to two hours of execution on mine machine and the author's experiments were executed with one hour of the time limit, then I can run my experiments with a time limit of two hours and have more precise results to compare.</p> <p>Hence, I would like to know if there is any way to calculate the equivalent time from two different machines (even if this time be approximated).</p> https://or.stackexchange.com/q/4238 3 Why does the getTime() function from cplex concert returns wrong value? Matheus Diógenes Andrade https://or.stackexchange.com/users/2862 2020-05-21T19:57:43Z 2021-11-20T14:58:34Z <p>I am running a MILP formulation (implemented in C++) with the Cplex Concert Technology 12.10, and I am trying to get the total elapsed time. So till the moment, I have tried three approaches: Be <code>cplex</code> my IloCplex object.</p> <ol> <li>Using the <code>clock()</code> function from the <code>time.h</code> C++ library:</li> </ol> <pre><code> cplex.setParam(IloCplex::Param::TimeLimit, time_limit); time_t start = clock(); cplex.solve(); double total_time = (double) (clock() - start) / (double) CLOCKS_PER_SEC; </code></pre> <ol start="2"> <li>Using the Cplex Concert function <code>getTime()</code>:</li> </ol> <pre><code> cplex.setParam(IloCplex::Param::TimeLimit, time_limit); cplex.solve(); double total_time = cplex.getTime(); </code></pre> <ol start="3"> <li>Using the Cplex Concert parameter <code>ClockType</code>:</li> </ol> <pre><code> cplex.setParam(IloCplex::Param::TimeLimit, time_limit); cplex.setParam(IloCplex::Param::ClockType, 2); cplex.solve(); double total_time = cplex.getTime(); </code></pre> <p>My code <strong>does not</strong> make use of any callbacks (informational callback neither), and does not set any value to the Cplex Concert parameter <code>Threads</code>, i.e., the piece of code <code>cplex.setParam(IloCplex::Param::Threads, n);</code> is not executed. Since such piece of code is not executed, then (according to the Cplex Concert official <a href="https://www.ibm.com/support/knowledgecenter/SSSA5P_12.8.0/ilog.odms.cplex.help/CPLEX/Parameters/topics/Threads.html" rel="nofollow noreferrer">documentation</a>) my code is using all available threads:</p> <blockquote> <p>When this parameter is at its default setting 0 (zero), and your application includes no callbacks or only an informational callback, CPLEX can use all available threads.</p> </blockquote> <p>The problem that I am facing is that in all these three approaches the <code>total_time</code> presents a wrong value. For example, if I set the variable <code>time_limit</code> to <code>120</code>, i.e., set the optimization time limit to <code>120</code> seconds, then the variable <code>total_time</code> presents a value much bigger than the expected one, such as <code>892.322</code>.</p> <p>So, I'd like to know if anyone already faced this problem before.</p>