0
$\begingroup$

Complete revision:

I have all 4 machines that can run in positive and negative directions which results in 2 outputs: $ P_1 = \begin{cases} -390 \le P_1 \le -300 & \text{for neg mode}\\ 0 & \text{for off mode}\\ 190 \le P_1 \le 380 & \text{for pos mode} \end{cases}$

These must match exactly the "Setpoint" while minimizing the efficiency. The efficiency for each machine is given by the equation: $ \frac{0.03}{190}x+0.937-0.2 \cdot i\\ \text{for i in the range number of machine} $ So far I tried formulating like this using 3 variables for each mode, but I still can't reformulate to constraint the output.

# Define the problem as a minimization problem, for each setpoint
for P_co in P_cos:
prob = LpProblem("HydroOptimization", LpMaximize)

# Define decision variables
P = [LpVariable(f"P{i}", P_min_p, P_max_t, cat="Continuous") for i in range(N)]
X = [LpVariable(f"X{i}{j}", cat ="Binary") for i in range(N) for j in ["p","t","z"]]
print(X)

Eff = LpVariable(f"Eff", lowBound=0, upBound= 1, cat='Continuous')



for i in range(N):
    prob += X[3*i] + X[3*i+1] + X[3*i+2] == 1, f"Pump Turbine Zero condition{i}"

    Less_eq = LpConstraint(e = P[i],
                            sense=LpConstraintLE,
                            rhs = P_max_p if X[3*i] == 1 else 0) 
    
    Great_eq = LpConstraint(e =P[i],
                            sense=LpConstraintGE,
                            rhs = P_min_t if X[3*i+1] == 1 else 0)
    
    Great_eq = LpConstraint(e =P[i],
                            sense=LpConstraintEQ,
                            rhs = 0 if X[3*i+2] == 1 else 0)
    
    prob.addConstraint(name = f"less_eq{i}", constraint=X[3*i]<=P[i]<=) 
    prob.addConstraint(name = f"Great_eq{i}", constraint=Great_eq) 

    prob.addConstraint(name = f"Prod{i}", constraint=  LpAffineExpression(P[i], X[3*i]) + LpAffineExpression(P[i], X[3*i+1]) + X[3*i+2]*0 == P[i])

print(prob)


prob += lpSum(P) == P_co, "Setpoint"

# Define objective function
for i in range(N):
    Eff += (-0.0000012*(P[i]-390)**2+1)*(P[i])
Eff = Eff/4*P_co
prob += Eff
print(prob)
# Solve the problem
status = prob.solve()

# Print the solution status
print("Solution status:", LpStatus[status])

# Print the decision variable values

print("P values:", [P[i].value() for i in range(N)])
print("Efficiency:",Eff.value())

Also later I would like to implement the start up and shutdown cost of each machine

Improve this question
$\endgroup$
5
  • $\begingroup$ @Reinderien the certain output is the setpoint, its an equality $\endgroup$
    – sino
    May 11 at 15:55
  • $\begingroup$ Between 0 and 600 usually around 400 $\endgroup$
    – sino
    May 11 at 21:46
  • $\begingroup$ It would be helpful if you can elaborate more on the question details, but it seems you want to set some conditional relations between efficiency and machines(?) If so, you need to imply the condition as $(x_1 \land \lnot x_2 \land \lnot x_3 \land \lnot x_4) \implies efficiency =0.8$ for the second row of your table. Is this what you are looking for? $\endgroup$
    – A.Omidi
    May 13 at 7:17
  • $\begingroup$ Yes exactly, sorry for the not so clear explanations, I'll try to edit it $\endgroup$
    – sino
    May 13 at 21:04
  • $\begingroup$ @sino I assume that you mean "maximize". I demonstrate a new answer with three variables per machine. $\endgroup$
    – Reinderien
    May 17 at 13:02

3 Answers 3

Reset to default
1
$\begingroup$

I assumed that the mode indicator is a variable otherwise, what you are looking for can be applied by a simple if-then loop on the constraints. To define the first row of your logical expression, it is equivalent to:

$$ Iff \quad (mod=-1) \implies (-390 \leq P_{1} \leq -300) $$

To linearize this constraint you need to define the binary auxiliary variables $z_{1}, \cdots, z_{3}$, and derive the following constraints:

$$ (mod = -1) \implies (z_{1} = 1) \iff (z_{1} = 0) \implies (mod \neq -1) $$ $$\cdots $$ $$ ((z_{2} = 1) \implies (-390 \leq P_{1})) \land ((z_{3} = 1) \implies (P_{1} \leq -300)) $$ $$\cdots $$ $$ z_{1} \implies (z_{2} \land z_{3}) $$ $$ \lnot z_{1} \lor (z_{2} \land z_{3}) $$ $$ (\lnot z_{1} \lor z_{2}) \land (\lnot z_{1} \lor z_{3}) $$ $$\cdots $$ $$ (1-z_{1}) + z_{2} \geq 1 \quad (1)$$ $$ (1-z_{1}) + z_{3} \geq 1 \quad (2)$$ $$ mod \geq 0 - M_1(z_{1}) \quad (3) $$ $$ mod \leq -2 + M_2(z_{1}) \quad (4) $$ $$ P_{1} \geq -390 - M_3(1-z_{2}) \quad (5) $$ $$ P_{1} \leq -300 + M_4(1-z_{3}) \quad (6) $$

Also, the rest parts can already be linearized in the above form.

Improve this answer
$\endgroup$
0
$\begingroup$

For the first 2 machines $m$ define a set of binary vars $p_{m,0}, p_{m,50},p_{m,100}$. Now for production variable $x \in R$ introduce constraints
(1) $ x_m = 50p_{m,50} + 100p_{m,100}$
$ (2) p_{m,0}+p_{m,50}+p_{m,100}=1 \quad \forall m =\{1,2\}$

For efficiency
$e_m = p_{m,50}-p_{m,100} $:
$e_m$ need not be an optimization variable. It could an expression to be used in the objective.

In fact if you want to reduce number of variables/constraints then for the 2 fixed machines $50p_{m,50}+100p_{m,100} $ can be used as production variable, instead of $x$ along with constraint (2) for each of fixed machines.

For the other two machines
(3) $b_m \le x_m \le 200b_m $
Here $b_m$ is a binary variable & can act as efficiency variable.

EDIT:
Based on new piecewise production function, you'd need the following constraints
$ -390\delta_1+190\delta_2 \le x \le -300\delta_1+380\delta_2$
$ \delta_1 + \delta_2 + \delta_3 = 1$: ensures at least one of the binaries is 1
where $ \delta \in \{0,1 \}^{3}$

Improve this answer
$\endgroup$
1
  • $\begingroup$ I also thought of using binary variables , but missed the condition of the sum =1. Will try to implement this on python $\endgroup$
    – sino
    May 12 at 12:41
0
$\begingroup$

Each machine needs at least two integral variables. I demonstrate here with one 'negative' variable and one 'positive' variable. I assume that in your description "minimize efficiency" is incorrect and you mean "maximize efficiency".

I also assume that the "set point" is the sum of the machine outputs, regardless of whether they're negative or positive. Your description left this ambiguous.

import pandas as pd
import pulp


def make_vars(index: int) -> pd.Series:
    output = pulp.LpVariable(
        name=f'P{index}', cat=pulp.LpContinuous)
    neg = pulp.LpVariable(name=f'neg{index}', cat=pulp.LpBinary)
    pos = pulp.LpVariable(name=f'pos{index}', cat=pulp.LpBinary)
    eff = 0.03/190*output + 0.937 - 0.2*index

    prob.addConstraint(name=f'excl{index}', constraint=neg + pos <= 1)

    '''
    sign  lo   hi
    -1  -390 -300
     0     0    0
    +1   190  380
    '''
    prob.addConstraint(name=f'lo{index}', constraint=output >= neg*-390 + pos*190)
    prob.addConstraint(name=f'hi{index}', constraint=output <= neg*-300 + pos*380)

    return pd.Series([output, neg, pos, eff], index=['output', 'neg', 'pos', 'eff'])


prob = pulp.LpProblem(name='HydroOptimization', sense=pulp.LpMaximize)
machines = pd.Series(
    pd.RangeIndex(start=0, stop=4), name='index'
).apply(make_vars)

set_point = 120
prob.addConstraint(
    name='set_point',
    constraint=pulp.lpSum(machines.output) == set_point)
prob.objective = pulp.lpSum(machines.eff)

print(machines)
print(prob)
prob.solve()
assert prob.status == pulp.LpStatusOptimal

machines[['output', 'neg', 'pos']] = machines[[
    'output', 'neg', 'pos'
]].applymap(pulp.LpVariable.value)
machines['eff'] = machines['eff'].apply(pulp.LpAffineExpression.value)
print(machines)

  output   neg   pos                           eff
0     P0  neg0  pos0  {P0: 0.00015789473684210527}
1     P1  neg1  pos1  {P1: 0.00015789473684210527}
2     P2  neg2  pos2  {P2: 0.00015789473684210527}
3     P3  neg3  pos3  {P3: 0.00015789473684210527}
HydroOptimization:
MAXIMIZE
0.00015789473684210527*P0 + 0.00015789473684210527*P1 + 0.00015789473684210527*P2 + 0.00015789473684210527*P3 + 2.548
SUBJECT TO
excl0: neg0 + pos0 <= 1

lo0: P0 + 390 neg0 - 190 pos0 >= 0

hi0: P0 + 300 neg0 - 380 pos0 <= 0

excl1: neg1 + pos1 <= 1

lo1: P1 + 390 neg1 - 190 pos1 >= 0

hi1: P1 + 300 neg1 - 380 pos1 <= 0

excl2: neg2 + pos2 <= 1

lo2: P2 + 390 neg2 - 190 pos2 >= 0

hi2: P2 + 300 neg2 - 380 pos2 <= 0

excl3: neg3 + pos3 <= 1

lo3: P3 + 390 neg3 - 190 pos3 >= 0

hi3: P3 + 300 neg3 - 380 pos3 <= 0

set_point: P0 + P1 + P2 + P3 = 120

VARIABLES
P0 free Continuous
P1 free Continuous
P2 free Continuous
P3 free Continuous
0 <= neg0 <= 1 Integer
0 <= neg1 <= 1 Integer
0 <= neg2 <= 1 Integer
0 <= neg3 <= 1 Integer
0 <= pos0 <= 1 Integer
0 <= pos1 <= 1 Integer
0 <= pos2 <= 1 Integer
0 <= pos3 <= 1 Integer

Welcome to the CBC MILP Solver 
Version: 2.10.3 
Build Date: Dec 15 2019 

At line 2 NAME          MODEL
At line 3 ROWS
At line 18 COLUMNS
At line 75 RHS
At line 89 BOUNDS
At line 102 ENDATA
Problem MODEL has 13 rows, 12 columns and 36 elements
Coin0008I MODEL read with 0 errors
Option for timeMode changed from cpu to elapsed
Continuous objective value is 0.0189474 - 0.00 seconds
Cgl0004I processed model has 12 rows, 11 columns (8 integer (8 of which binary)) and 36 elements
Cbc0038I Initial state - 1 integers unsatisfied sum - 0.315789
Cbc0038I Pass   1: suminf.    0.31579 (1) obj. -0.0189474 iterations 1
Cbc0038I Pass   2: suminf.    0.17949 (1) obj. -0.0189474 iterations 5
Cbc0038I Pass   3: suminf.    0.10526 (1) obj. -0.0189474 iterations 6
Cbc0038I Pass   4: suminf.    0.00000 (0) obj. -0.0189474 iterations 4
Cbc0038I Solution found of -0.0189474
Cbc0038I Relaxing continuous gives -0.0189474
Cbc0038I Before mini branch and bound, 4 integers at bound fixed and 0 continuous
Cbc0038I Mini branch and bound did not improve solution (0.00 seconds)
Cbc0038I After 0.00 seconds - Feasibility pump exiting with objective of -0.0189474 - took 0.00 seconds
Cbc0012I Integer solution of -0.018947368 found by feasibility pump after 0 iterations and 0 nodes (0.00 seconds)
Cbc0001I Search completed - best objective -0.018947368421052, took 0 iterations and 0 nodes (0.00 seconds)
Cbc0035I Maximum depth 0, 0 variables fixed on reduced cost
Cuts at root node changed objective from -0.0189474 to -0.0189474
Probing was tried 0 times and created 0 cuts of which 0 were active after adding rounds of cuts (0.000 seconds)
Gomory was tried 0 times and created 0 cuts of which 0 were active after adding rounds of cuts (0.000 seconds)
Knapsack was tried 0 times and created 0 cuts of which 0 were active after adding rounds of cuts (0.000 seconds)
Clique was tried 0 times and created 0 cuts of which 0 were active after adding rounds of cuts (0.000 seconds)
MixedIntegerRounding2 was tried 0 times and created 0 cuts of which 0 were active after adding rounds of cuts (0.000 seconds)
FlowCover was tried 0 times and created 0 cuts of which 0 were active after adding rounds of cuts (0.000 seconds)
TwoMirCuts was tried 0 times and created 0 cuts of which 0 were active after adding rounds of cuts (0.000 seconds)
ZeroHalf was tried 0 times and created 0 cuts of which 0 were active after adding rounds of cuts (0.000 seconds)

Result - Optimal solution found

Objective value:                0.01894737
Enumerated nodes:               0
Total iterations:               0
Time (CPU seconds):             0.00
Time (Wallclock seconds):       0.00

Option for printingOptions changed from normal to all
Total time (CPU seconds):       0.00   (Wallclock seconds):       0.01

   output  neg  pos       eff
0  -300.0  1.0  0.0  0.889632
1   190.0  0.0  1.0  0.767000
2     0.0  0.0  0.0  0.537000
3   230.0  0.0  1.0  0.373316
Improve this answer
$\endgroup$
2
  • $\begingroup$ Would you also use binary coeffiencts if the fixed machine were to run at values that arent multiples in between like [0,70,123]? $\endgroup$
    – sino
    May 12 at 13:13
  • $\begingroup$ Yes, you would use binary variables, but the setup of the problem would change because you'd no longer be able to use a single integer and scale it. The size of the problem would increase. $\endgroup$
    – Reinderien
    May 12 at 15:06

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.