Skip to content
Snippets Groups Projects
Commit f2ecd0b0 authored by Steinmann's avatar Steinmann
Browse files

started implemented Flow constraint

parent 71c78a82
Branches
No related tags found
No related merge requests found
Show whitespace changes
Inline Side-by-side
optimierung_pymoo.ipynb
+ 54
− 100
View file @ f2ecd0b0
%% Cell type:markdown id: tags: %% Cell type:markdown id: tags:
Formulieren der Optimierungsgleichung in pymoo Formulieren der Optimierungsgleichung in pymoo
%% Cell type:markdown id: tags: %% Cell type:markdown id: tags:
Es gilt die Kontinuitätsgleichung: Es gilt die Kontinuitätsgleichung:
$ \Sigma \dot{V}_k(t) = O$ $ \Sigma \dot{V}_k(t) = O$
und die aus der Topologie resultierende Inzidenzmatrix $A_i$ und die aus der Topologie resultierende Inzidenzmatrix $A_i$
sowie die aus dem Pumpenkennfeld folgende Beziehung: sowie die aus dem Pumpenkennfeld folgende Beziehung:
$\Delta p=\alpha_1 Q^2+\alpha_2 Q n+\alpha_3 n^2 : n \in \{0\} \cup [n_{\mathrm{min}},n_{\mathrm{max}}] $ $\Delta p=\alpha_1 Q^2+\alpha_2 Q n+\alpha_3 n^2 : n \in \{0\} \cup [n_{\mathrm{min}},n_{\mathrm{max}}] $
$P=\beta_1 Q^3+\beta_2 Q^2 n+\beta_3 Q n^2+\beta_4n^3+\beta_5$ $P=\beta_1 Q^3+\beta_2 Q^2 n+\beta_3 Q n^2+\beta_4n^3+\beta_5$
und die beziehung für den Druckverlust an den Ventilen: und die beziehung für den Druckverlust an den Ventilen:
$\Delta p_{\mathrm{loss}} = - \frac{1}{2} \varrho \zeta \left(\frac{Q}{A}\right)^2 = -l Q^2 :l\in [l_{\mathrm{min}}:\infty )$ $\Delta p_{\mathrm{loss}} = - \frac{1}{2} \varrho \zeta \left(\frac{Q}{A}\right)^2 = -l Q^2 :l\in [l_{\mathrm{min}}:\infty )$
nun soll für einen Gegebenen Volumenstrom $Q$ eine Optimale Drehzahl bestimmt werden, welche die Pumpenlesitung minimiert. nun soll für einen Gegebenen Volumenstrom $Q$ eine Optimale Drehzahl bestimmt werden, welche die Pumpenlesitung minimiert.
$$ $$
\begin{align*} \begin{align*}
\mathrm{min} \sum_{p \in \mathcal{P}} Po_{p} \\ \mathrm{min} \sum_{p \in \mathcal{P}} Po_{p} \\
Q_p = Q_{\mathrm{target}, i} \\ Q_{p,i} \geq \sum_{strang} Q_v + \sum_{strang} Q_p \\
Q_p Q_p , n\epsilon [n_{min},n_{max}] \\
, n\epsilon [n_{min},n_{max}] \\
\overrightarrow{n} = (1,n,n^2,n^3)^T \\ \overrightarrow{n} = (1,n,n^2,n^3)^T \\
min P = A \overrightarrow{n} \\ min P = A \overrightarrow{n} \\
-n\leq n_{min} \\ -n\leq n_{min} \\
n\leq n_{max} n\leq n_{max}
\end{align*} \end{align*}
$$ $$
Förderhöhe als constraint continuität fomulieren pro strang Förderhöhe als constraint continuität fomulieren pro strang
%% Cell type:code id: tags: %% Cell type:code id: tags:
``` python ``` python
!pip install pyomo !pip install pyomo
``` ```
%% Cell type:code id: tags: %% Cell type:code id: tags:
``` python ``` python
#Pump-Powercurve and Pump-Hightcurve #Pump-Powercurve and Pump-Hightcurve
import regression_own import regression_own
(LR_H,LR_P)=regression_own.regress_pump() (LR_H,LR_P)=regression_own.regress_pump()
``` ```
%% Output %% Output
R^20.9998289611292903 R^20.9998289611292903
R^20.9994449560888792 R^20.9994449560888792
%% Cell type:code id: tags: %% Cell type:code id: tags:
``` python ``` python
#Graph constroctor #Graph constroctor
#Alle Ventile sind direkt mit der Quelle/Senke Verbunden #Alle Ventile sind direkt mit der Quelle/Senke Verbunden
import multiDiGraph as gr import multiDiGraph as gr
nodes =['source','pump1','pump2','valveA','valveB','valveC'] nodes =['source','pump1','pump2','valveA','valveB','valveC']
graph = gr.construct_graph('source',('source','pump1',0.),('pump1','pump2',0.),('pump2','valveA',0.),('pump2','valveB',0.), graph = gr.construct_graph('source',('source','pump1',0.),('pump1','pump2',0.),('pump2','valveA',0.),('pump2','valveB',0.),
('pump1','valveC',0.),('valveA','source',4.),('valveB','source',4.),('valveC','source',4.)) ('pump1','valveC',0.),('valveA','source',4.),('valveB','source',4.),('valveC','source',4.))
#ist das notwendig?!? #ist das notwendig?!?
for node in graph.nodes: for node in graph.nodes:
#definieren der Drehzahl für jede Pumpe im graphen #definieren der Drehzahl für jede Pumpe im graphen
#inizieren des Durchflusses für jedes Ventil im Graphen #inizieren des Durchflusses für jedes Ventil im Graphen
if 'pump' in node: if 'pump' in node:
graph.nodes[node]['n']=750/3600 graph.nodes[node]['n']=750/3600
else: else:
graph.nodes[node]['n']=None graph.nodes[node]['n']=None
graph.nodes[node]['flow']=0. graph.nodes[node]['flow']=0.
if 'valve' in node: if 'valve' in node:
graph.nodes[node]['flow']= graph[node]['source'][0]['weight'] graph.nodes[node]['flow']= graph[node]['source'][0]['weight']
for node in graph.nodes: for node in graph.nodes:
#Berechnen des Durchflusses im Knoten #Berechnen des Durchflusses im Knoten
if 'valve' in node: if 'valve' in node:
continue continue
for inF in graph.predecessors(node): for inF in graph.predecessors(node):
graph.nodes[node]['flow'] += graph[inF][node][0]['weight'] graph.nodes[node]['flow'] += graph[inF][node][0]['weight']
#Berechnen des Durchflusses der abgehenden Kanten #Berechnen des Durchflusses der abgehenden Kanten
tempF=graph.nodes[node]['flow'] tempF=graph.nodes[node]['flow']
SC=0 SC=0
for outF in graph.successors(node): for outF in graph.successors(node):
if 'valve' in outF: if 'valve' in outF:
graph[node][outF][0]['weight']=graph.nodes[outF]['flow'] graph[node][outF][0]['weight']=graph.nodes[outF]['flow']
tempF=tempF - graph.nodes[outF]['flow'] tempF=tempF - graph.nodes[outF]['flow']
else: else:
SC+=1 SC+=1
for outF in graph.successors(node): for outF in graph.successors(node):
if SC!=0. and not'valve' in outF: if SC!=0. and not'valve' in outF:
graph[node][outF][0]['weight']=tempF/SC graph[node][outF][0]['weight']=tempF/SC
else:continue else:continue
print(graph.nodes.data('flow')) print(graph.nodes.data('flow'))
``` ```
%% Output %% Output
[('source', 12.0), ('pump1', 12.0), ('pump2', 8.0), ('valveA', 4.0), ('valveB', 4.0), ('valveC', 4.0)] [('source', 12.0), ('pump1', 12.0), ('pump2', 8.0), ('valveA', 4.0), ('valveB', 4.0), ('valveC', 4.0)]
%% Cell type:code id: tags: %% Cell type:code id: tags:
``` python ``` python
import networkx as nx import networkx as nx
Mtrx= nx.incidence_matrix(graph,nodes,oriented=True) Mtrx= nx.incidence_matrix(graph,nodes,oriented=True)
``` ```
%% Output
<Compressed Sparse Column sparse array of dtype 'float64'
with 16 stored elements and shape (6, 8)>
Coords Values
(0, 0) -1.0
(1, 0) 1.0
(1, 1) -1.0
(2, 1) 1.0
(1, 2) -1.0
(5, 2) 1.0
(2, 3) -1.0
(3, 3) 1.0
(2, 4) -1.0
(4, 4) 1.0
(0, 5) 1.0
(3, 5) -1.0
(0, 6) 1.0
(4, 6) -1.0
(0, 7) 1.0
(5, 7) -1.0
%% Cell type:code id: tags: %% Cell type:code id: tags:
``` python ``` python
import networkx as nx import networkx as nx
def create_dict(GR:nx.multidigraph): def create_dict(GR:nx.multidigraph):
data={None:{'nodes':{}, data={None:{'nodes':{},
'pumps':{}, 'pumps':{},
'valves':{}, 'valves':{},
} }
} }
for node in GR.nodes: for node in GR.nodes:
data[None]['nodes'][node]=None data[None]['nodes'][node]=None
data[None]['Q'][node]=GR.nodes[node]['flow'] data[None]['Q'][node]=GR.nodes[node]['flow']
if 'pump' in node: if 'pump' in node:
data[None]['pumps'][node]=None data[None]['pumps'][node]=None
data[None]['n'][node]=0. data[None]['n'][node]=0.
if 'valve' in node: if 'valve' in node:
data[None]['valves'][node]=None data[None]['valves'][node]=None
return data return data
``` ```
%% Cell type:markdown id: tags: %% Cell type:markdown id: tags:
Durchfluss aus Incidenzmatrix beerechnen Durchfluss aus Incidenzmatrix beerechnen
Zeilen = knoten Zeilen = knoten
Spalten = kanten Spalten = kanten
Summe pro knoten = 0 Summe pro knoten = 0
Q pump muss größer gleich sein als alle nachfolgenden durchflüsse
%% Cell type:code id: tags: %% Cell type:code id: tags:
``` python ``` python
#defining abstract modell for given Network #defining abstract modell for given Network
import pyomo.environ as pyo import pyomo.environ as pyo
from pyomo.dataportal import DataPortal from pyomo.dataportal import DataPortal
import numpy as np import numpy as np
from sklearn.linear_model import LinearRegression from sklearn.linear_model import LinearRegression
modell = pyo.AbstractModel() modell = pyo.AbstractModel()
#notwendige Mengen zur Berechnung der Constraints #notwendige Mengen zur Berechnung der Constraints
modell.nodes = pyo.Set() modell.nodes = pyo.Set()
modell.pumps = pyo.Set() modell.pumps = pyo.Set()
modell.valves = pyo.Set() modell.valves = pyo.Set()
modell.Q_valve=pyo.Param(modell.valves)
#Optimierungsvariable #Optimierungsvariable
modell.n = pyo.Var(modell.pumps,bounds=(0,1)) modell.n = pyo.Var(modell.pumps,bounds=(750/3600,1))
modell.Q = pyo.Var(modell.nodes) modell.Q_pump = pyo.Var(modell.pumps)
#Objective #Objective
def PumpPower(modell): def PumpPower(modell):
return sum(np.dot( return sum(np.dot(
np.array( np.array(
[modell.Q[i]**3,(modell.Q[i]**2)*modell.n[i],modell.Q[i]*modell.n[i]**2,modell.n[i]**3] [modell.Q_pump[i]**3,(modell.Q_pump[i]**2)*modell.n[i],modell.Q_pump[i]*modell.n[i]**2,modell.n[i]**3]
),LR_P.coef_ ),LR_P.coef_
) for i in modell.pumps) ) for i in modell.pumps)
modell.Power_Objective = pyo.Objective(rule=PumpPower,sense=pyo.minimize) modell.Power_Objective = pyo.Objective(rule=PumpPower,sense=pyo.minimize)
def PumpFlow(modell,pump):
pump=np.dot(np.array([modell.Q_pump[pump]**2,modell.n[pump]*modell.Q_pump[pump],modell.n[pump]**2]),LR_H.coef_)
return pump>=sum(modell.Q_valve[node] for node in graph.successors(pump) if node in modell.valves)+sum(modell.Q_pump[n] for n in graph.successors(node) if node in modell.pumps)
modell.Flow_Objective = pyo.Objective(rule=PumpFlow,sense=pyo.as_boolean)
#Constaints #Constaints
def continuityRule(modell,node): def continuityRule(modell,node):
return sum(modell.Q[i] for i in graph.predecessors(node))==sum(modell.Q[j] for j in graph.successors(node)) return sum(modell.Q_pump[i] for i in graph.predecessors(node))==sum(modell.Q_pump[j] for j in graph.successors(node))
#alternative #alternative
def continuityRule2(modell,node): def continuityRule2(modell,node):
return 0.==sum(graph[node][i][0]['weight'] for i in graph[node]) return 0.==sum(graph[node][i][0]['weight'] for i in graph[node])
#continuity adjustment for change in hight needed #continuity adjustment for change in hight needed
#construction of test Data dictionairy missing #construction of test Data dictionairy missing
TestData={ TestData={
None:{ None:{
'nodes':[key for key in graph.nodes.keys()], 'Q_valve':{'valveA':4.,'valveB':4.,'valveC':4.},
'pumps':[key for key in graph.nodes.keys() if 'pump' in key], 'nodes':[key for key in graph.nodes.keys()],
'valves':[key for key in graph.nodes.keys() if 'valve' in key], 'pumps':[key for key in graph.nodes.keys() if 'pump' in key],
'valves':[key for key in graph.nodes.keys() if 'valve' in key],
} }
} }
print(TestData) print(TestData)
#data=DataPortal(data_dict=TestData) #data=DataPortal(data_dict=TestData)
#Optimierungsgleichung #Optimierungsgleichung
#modell.pump_constraint = pyo.Constraint(expr=sum(modell.nodes[k] for k in modell.nodes)==0,rule=continuityRule) #modell.pump_constraint = pyo.Constraint(expr=sum(modell.nodes[k] for k in modell.nodes)==0,rule=continuityRule)
#instance=modell.create_instance(graph,LR_H) #instance=modell.create_instance(graph,LR_H)
#instance.obj = pyo.Objective(expr=sum(PumpPower(modell.Q[i],modell.n[i],LR_P) for i in modell.pumps),sense=min) #instance.obj = pyo.Objective(expr=sum(PumpPower(modell.Q_pump[i],modell.n[i],LR_P) for i in modell.pumps),sense=min)
``` ```
%% Output %% Output
{None: {'nodes': ['source', 'pump1', 'pump2', 'valveA', 'valveB', 'valveC'], 'pumps': ['pump1', 'pump2'], 'valves': ['valveA', 'valveB', 'valveC']}} {None: {'Q_valve': {'valveA': 4.0, 'valveB': 4.0, 'valveC': 4.0}, 'nodes': ['source', 'pump1', 'pump2', 'valveA', 'valveB', 'valveC'], 'pumps': ['pump1', 'pump2'], 'valves': ['valveA', 'valveB', 'valveC']}}
%% Cell type:markdown id: tags: %% Cell type:markdown id: tags:
Frage: gibt es nur eine Lösung für Drehzahl? Frage: gibt es nur eine Lösung für Drehzahl?
Bsp. Optimierung nach Dezentraler Pumpe um modell zu prüfen Bsp. Optimierung nach Dezentraler Pumpe um modell zu prüfen
%% Cell type:code id: tags: %% Cell type:code id: tags:
``` python ``` python
from pyomo.opt import SolverFactory from pyomo.opt import SolverFactory
opt = pyo.SolverFactory('scipampl', executable=r'C:\Program Files\SCIPOptSuite 9.2.0\bin\scip.exe') opt = pyo.SolverFactory('scipampl', executable=r'C:\Program Files\SCIPOptSuite 9.2.0\bin\scip.exe')
instance = modell.create_instance(TestData) instance = modell.create_instance(TestData)
result=opt.solve(instance, tee=True) result=opt.solve(instance, tee=True)
``` ```
%% Output %% Output
SCIP version 9.2.0 [precision: 8 byte] [memory: block] [mode: optimized] [LP solver: Soplex 7.1.2] [GitHash: 74cea9222e] ERROR: Rule failed when generating expression for Objective Flow_Objective
Copyright (c) 2002-2024 Zuse Institute Berlin (ZIB) with index None: KeyError: "Index 'None' is not valid for indexed component
'Q_pump'"
External libraries: ERROR: Constructing component 'Flow_Objective' from data=None failed:
Soplex 7.1.2 Linear Programming Solver developed at Zuse Institute Berlin (soplex.zib.de) [GitHash: b040369c] KeyError: "Index 'None' is not valid for indexed component 'Q_pump'"
CppAD 20180000.0 Algorithmic Differentiation of C++ algorithms developed by B. Bell (github.com/coin-or/CppAD)
TinyCThread 1.2 small portable implementation of the C11 threads API (tinycthread.github.io) ---------------------------------------------------------------------------
MPIR 3.0.0 Multiple Precision Integers and Rationals Library developed by W. Hart (mpir.org) KeyError Traceback (most recent call last)
ZIMPL 3.6.2 Zuse Institute Mathematical Programming Language developed by T. Koch (zimpl.zib.de) Cell In[19], line 4
AMPL/MP 690e9e7 AMPL .nl file reader library (github.com/ampl/mp) 1 from pyomo.opt import SolverFactory
PaPILO 2.4.0 parallel presolve for integer and linear optimization (github.com/scipopt/papilo) (built with TBB) [GitHash: 2d9fe29f] 3 opt = pyo.SolverFactory('scipampl', executable=r'C:\Program Files\SCIPOptSuite 9.2.0\bin\scip.exe')
Nauty 2.8.8 Computing Graph Automorphism Groups by Brendan D. McKay (users.cecs.anu.edu.au/~bdm/nauty) ----> 4 instance = modell.create_instance(TestData)
sassy 1.1 Symmetry preprocessor by Markus Anders (github.com/markusa4/sassy) 6 result=opt.solve(instance, tee=True)
Ipopt 3.14.16 Interior Point Optimizer developed by A. Waechter et.al. (github.com/coin-or/Ipopt) File ~\AppData\Roaming\Python\Python312\site-packages\pyomo\core\base\PyomoModel.py:734, in Model.create_instance(self, filename, data, name, namespace, namespaces, profile_memory, report_timing, **kwds)
731 if None not in _namespaces:
user parameter file <scip.set> not found - using default parameters 732 _namespaces.append(None)
read problem <C:\Users\STEINM~1\AppData\Local\Temp\tmptn91lkdj.pyomo.nl> --> 734 instance.load(data, namespaces=_namespaces, profile_memory=profile_memory)
============ 736 #
737 # Indicate that the model is concrete/constructed
original problem has 5 variables (0 bin, 0 int, 0 impl, 5 cont) and 1 constraints 738 #
739 instance._constructed = True
solve problem File ~\AppData\Roaming\Python\Python312\site-packages\pyomo\core\base\PyomoModel.py:771, in Model.load(self, arg, namespaces, profile_memory)
============= 769 msg = "Cannot load model model data from with object of type '%s'"
770 raise ValueError(msg % str(type(arg)))
feasible solution found by trivial heuristic after 0.0 seconds, objective value 0.000000e+00 --> 771 self._load_model_data(dp, namespaces, profile_memory=profile_memory)
presolving: File ~\AppData\Roaming\Python\Python312\site-packages\pyomo\core\base\PyomoModel.py:823, in Model._load_model_data(self, modeldata, namespaces, **kwds)
(0.0s) symmetry computation started: requiring (bin +, int +, cont +), (fixed: bin -, int -, cont -) 820 if component.ctype is Model:
(0.0s) symmetry computation finished: 1 generators found (max: 1500, log10 of symmetry group size: 0.0) (symcode time: 0.00) 821 continue
dynamic symmetry handling statistics: --> 823 self._initialize_component(
orbitopal reduction: no components 824 modeldata, namespaces, component_name, profile_memory
orbital reduction: no components 825 )
lexicographic reduction: 1 permutations with support sizes 4 827 # Note: As is, connectors are expanded when using command-line pyomo but not calling model.create(...) in a Python script.
handled 1 out of 1 symmetry components 828 # John says this has to do with extension points which are called from commandline but not when writing scripts.
presolving (1 rounds: 1 fast, 1 medium, 1 exhaustive): 829 # Uncommenting the next two lines switches this (command-line fails because it tries to expand connectors twice)
0 deleted vars, 0 deleted constraints, 0 added constraints, 0 tightened bounds, 0 added holes, 0 changed sides, 0 changed coefficients 830 # connector_expander = ConnectorExpander()
0 implications, 0 cliques 831 # connector_expander.apply(instance=self)
presolved problem has 5 variables (0 bin, 0 int, 0 impl, 5 cont) and 1 constraints 833 if profile_memory >= 2 and pympler_available:
1 constraints of type <nonlinear> File ~\AppData\Roaming\Python\Python312\site-packages\pyomo\core\base\PyomoModel.py:871, in Model._initialize_component(self, modeldata, namespaces, component_name, profile_memory)
Presolving Time: 0.00 863 logger.debug(
transformed 3/3 original solutions to the transformed problem space 864 "Constructing %s '%s' on %s from data=%s",
865 declaration.__class__.__name__,
(...)
%% Cell type:code id: tags: 868 str(data),
869 )
``` python 870 try:
opt --> 871 declaration.construct(data)
``` 872 except:
873 err = sys.exc_info()[1]
%% Cell type:code id: tags: File ~\AppData\Roaming\Python\Python312\site-packages\pyomo\core\base\objective.py:335, in Objective.construct(self, data)
333 # Bypass the index validation and create the member directly
``` python 334 for index in self.index_set():
``` --> 335 ans = self._setitem_when_not_present(index, rule(block, index))
336 if ans is not None:
337 ans.set_sense(self._init_sense(block, index))
File ~\AppData\Roaming\Python\Python312\site-packages\pyomo\core\base\initializer.py:349, in IndexedCallInitializer.__call__(self, parent, idx)
347 return self._fcn(parent, *idx)
348 else:
--> 349 return self._fcn(parent, idx)
Cell In[18], line 30, in PumpFlow(modell, pump)
29 def PumpFlow(modell,pump):
---> 30 pump=np.dot(np.array([modell.Q_pump[pump]**2,modell.n[pump]*modell.Q_pump[pump],modell.n[pump]**2]),LR_H.coef_)
31 return pump>=sum(modell.Q_valve[node] for node in graph.successors(pump) if node in modell.valves)+sum(modell.Q_pump[n] for n in graph.successors(node) if node in modell.pumps)
File ~\AppData\Roaming\Python\Python312\site-packages\pyomo\core\base\var.py:999, in IndexedVar.__getitem__(self, args)
997 def __getitem__(self, args) -> VarData:
998 try:
--> 999 return super().__getitem__(args)
1000 except RuntimeError:
1001 tmp = args if args.__class__ is tuple else (args,)
File ~\AppData\Roaming\Python\Python312\site-packages\pyomo\core\base\indexed_component.py:648, in IndexedComponent.__getitem__(self, index)
646 if isinstance(index, EXPR.GetItemExpression):
647 return index
--> 648 validated_index = self._validate_index(index)
649 if validated_index is not index:
650 index = validated_index
File ~\AppData\Roaming\Python\Python312\site-packages\pyomo\core\base\indexed_component.py:870, in IndexedComponent._validate_index(self, idx)
863 raise KeyError(
864 "Cannot treat the scalar component '%s' "
865 "as an indexed component" % (self.name,)
866 )
867 #
868 # Raise an exception
869 #
--> 870 raise KeyError(
871 "Index '%s' is not valid for indexed component '%s'"
872 % (normalized_idx, self.name)
873 )
KeyError: "Index 'None' is not valid for indexed component 'Q_pump'"
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment