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

finished solver implementation

parent c1674030
No related branches found
No related tags found
No related merge requests found
Show whitespace changes
Inline Side-by-side
optimierung_pymoo.ipynb
+ 128
− 56
View file @ 71c78a82
%% Cell type:markdown id: tags:
Formulieren der Optimierungsgleichung in pymoo
%% Cell type:markdown id: tags:
Es gilt die Kontinuitätsgleichung:
$ \Sigma \dot{V}_k(t) = O$
und die aus der Topologie resultierende Inzidenzmatrix $A_i$
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}}] $
$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:
$\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.
$$
\begin{align*}
\mathrm{min} \sum_{p \in \mathcal{P}} Po_{p} \\
Q_p = Q_{\mathrm{target}, i} \\
Q_p
, n\epsilon [n_{min},n_{max}] \\
\overrightarrow{n} = (1,n,n^2,n^3)^T \\
min P = A \overrightarrow{n} \\
-n\leq n_{min} \\
n\leq n_{max}
\end{align*}
$$
Förderhöhe als constraint continuität fomulieren pro strang
%% Cell type:code id: tags:
``` python
!pip install pyomo
```
%% Cell type:code id: tags:
``` python
#Pump-Powercurve and Pump-Hightcurve
import regression_own
(LR_H,LR_P)=regression_own.regress_pump()
```
%% Output
R^20.9998289611292903
R^20.9994449560888792
%% Cell type:code id: tags:
``` python
#Graph constroctor
#Alle Ventile sind direkt mit der Quelle/Senke Verbunden
import multiDiGraph as gr
nodes =['source','pump1','pump2','valveA','valveB','valveC']
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.))
#ist das notwendig?!?
for node in graph.nodes:
#definieren der Drehzahl für jede Pumpe im graphen
#inizieren des Durchflusses für jedes Ventil im Graphen
if 'pump' in node:
graph.nodes[node]['n']=0.
graph.nodes[node]['n']=750/3600
else:
graph.nodes[node]['n']=None
graph.nodes[node]['flow']=0.
if 'valve' in node:
graph.nodes[node]['flow']= graph[node]['source'][0]['weight']
for node in graph.nodes:
#Berechnen des Durchflusses im Knoten
if 'valve' in node:
continue
for inF in graph.predecessors(node):
graph.nodes[node]['flow'] += graph[inF][node][0]['weight']
#Berechnen des Durchflusses der abgehenden Kanten
tempF=graph.nodes[node]['flow']
SC=0
for outF in graph.successors(node):
if 'valve' in outF:
graph[node][outF][0]['weight']=graph.nodes[outF]['flow']
tempF=tempF - graph.nodes[outF]['flow']
else:
SC+=1
for outF in graph.successors(node):
if SC!=0. and not'valve' in outF:
graph[node][outF][0]['weight']=tempF/SC
else:continue
print(graph.nodes.data('flow'))
```
%% Output
[('source', 12.0), ('pump1', 12.0), ('pump2', 8.0), ('valveA', 4.0), ('valveB', 4.0), ('valveC', 4.0)]
%% Cell type:code id: tags:
``` python
import networkx as nx
Mtrx= nx.incidence_matrix(graph,nodes,oriented=True)
print(Mtrx)
```
%% 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
(0, 1) 1.0
(3, 1) -1.0
(0, 2) 1.0
(4, 2) -1.0
(0, 3) 1.0
(5, 3) -1.0
(1, 4) -1.0
(2, 4) 1.0
(1, 5) -1.0
(5, 5) 1.0
(2, 6) -1.0
(3, 6) 1.0
(2, 7) -1.0
(4, 7) 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:
``` python
import networkx as nx
def create_dict(GR:nx.multidigraph):
data={None:{'nodes':{},
'pumps':{},
'valves':{},
}
}
for node in GR.nodes:
data[None]['nodes'][node]=None
data[None]['Q'][node]=GR.nodes[node]['flow']
if 'pump' in node:
data[None]['pumps'][node]=None
data[None]['n'][node]=0.
if 'valve' in node:
data[None]['valves'][node]=None
return data
```
%% Cell type:markdown id: tags:
Durchfluss aus Incidenzmatrix beerechnen
Zeilen = knoten
Spalten = kanten
Summe pro knoten = 0
%% Cell type:code id: tags:
``` python
#defining abstract modell for given Network
import pyomo.environ as pyo
from pyomo.dataportal import DataPortal
import numpy as np
from sklearn.linear_model import LinearRegression
modell = pyo.AbstractModel()
#notwendige Mengen zur Berechnung der Constraints
modell.nodes = pyo.Set()
modell.pumps = pyo.Set()
modell.valves = pyo.Set()
#Optimierungsvariable
modell.n = pyo.Var(modell.pumps,bounds=(750/3600,1))
modell.n = pyo.Var(modell.pumps,bounds=(0,1))
modell.Q = pyo.Var(modell.nodes)
#Objective
def PumpPower(modell):
return sum(np.dot(np.array([modell.Q[i]**3,(modell.Q[i]**2)*modell.n[i],modell.Q[i]*modell.n[i]**2,modell.n[i]**3]),LR_P.coef_) for i in modell.pumps)
modell.Power_Objective = pyo.Objective(rule=PumpPower,sense=min)
return sum(np.dot(
np.array(
[modell.Q[i]**3,(modell.Q[i]**2)*modell.n[i],modell.Q[i]*modell.n[i]**2,modell.n[i]**3]
),LR_P.coef_
) for i in modell.pumps)
modell.Power_Objective = pyo.Objective(rule=PumpPower,sense=pyo.minimize)
#Constaints
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))
#alternative
def continuityRule2(modell,node):
return 0.==sum(graph[node][i][0]['weight'] for i in graph[node])
#continuity adjustment for change in hight needed
#construction of test Data dictionairy missing
dict_Q={}
dict_n={}
for key in graph.nodes.keys():
dict_Q[key]=graph.nodes[key]['flow']
dict_n[key]=graph.nodes[key]['n']
TestData={
None:{
'nodes':[key for key in graph.nodes.keys()],
'pumps':[key for key in graph.nodes.keys() if 'pump' in key],
'valves':[key for key in graph.nodes.keys() if 'valve' in key],
'Q': dict_Q,
'n': dict_n
}
}
print(TestData)
data=DataPortal(data_dict=TestData)
instance = modell.create_instance(data)
#data=DataPortal(data_dict=TestData)
#Optimierungsgleichung
#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.obj = pyo.Objective(expr=sum(PumpPower(modell.Q[i],modell.n[i],LR_P) for i in modell.pumps),sense=min)
```
%% Output
{None: {'nodes': ['source', 'pump1', 'pump2', 'valveA', 'valveB', 'valveC'], 'pumps': ['pump1', 'pump2'], 'valves': ['valveA', 'valveB', 'valveC'], 'Q': {'source': 12.0, 'pump1': 12.0, 'pump2': 8.0, 'valveA': 8.0, 'valveB': 8.0, 'valveC': 8.0}, 'n': {'source': None, 'pump1': 0.0, 'pump2': 0.0, 'valveA': None, 'valveB': None, 'valveC': None}}}
ERROR: Rule failed when initializing variable for Var n with index None:
AssertionError:
ERROR: Constructing component 'n' from data={'source': None, 'pump1': 0.0,
'pump2': 0.0, 'valveA': None, 'valveB': None, 'valveC': None} failed:
AssertionError:
---------------------------------------------------------------------------
AssertionError Traceback (most recent call last)
Cell In[76], line 44
42 print(TestData)
43 data=DataPortal(data_dict=TestData)
---> 44 instance = modell.create_instance(data)
46 #Optimierungsgleichung
47 #modell.pump_constraint = pyo.Constraint(expr=sum(modell.nodes[k] for k in modell.nodes)==0,rule=continuityRule)
48 #instance=modell.create_instance(graph,LR_H)
49 #instance.obj = pyo.Objective(expr=sum(PumpPower(modell.Q[i],modell.n[i],LR_P) for i in modell.pumps),sense=min)
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:
732 _namespaces.append(None)
--> 734 instance.load(data, namespaces=_namespaces, profile_memory=profile_memory)
736 #
737 # Indicate that the model is concrete/constructed
738 #
739 instance._constructed = True
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)))
--> 771 self._load_model_data(dp, namespaces, profile_memory=profile_memory)
File ~\AppData\Roaming\Python\Python312\site-packages\pyomo\core\base\PyomoModel.py:823, in Model._load_model_data(self, modeldata, namespaces, **kwds)
820 if component.ctype is Model:
821 continue
--> 823 self._initialize_component(
824 modeldata, namespaces, component_name, profile_memory
825 )
827 # Note: As is, connectors are expanded when using command-line pyomo but not calling model.create(...) in a Python script.
828 # John says this has to do with extension points which are called from commandline but not when writing scripts.
829 # Uncommenting the next two lines switches this (command-line fails because it tries to expand connectors twice)
830 # connector_expander = ConnectorExpander()
831 # connector_expander.apply(instance=self)
833 if profile_memory >= 2 and pympler_available:
File ~\AppData\Roaming\Python\Python312\site-packages\pyomo\core\base\PyomoModel.py:871, in Model._initialize_component(self, modeldata, namespaces, component_name, profile_memory)
863 logger.debug(
864 "Constructing %s '%s' on %s from data=%s",
865 declaration.__class__.__name__,
(...)
868 str(data),
869 )
870 try:
--> 871 declaration.construct(data)
872 except:
873 err = sys.exc_info()[1]
File ~\AppData\Roaming\Python\Python312\site-packages\pyomo\core\base\var.py:735, in Var.construct(self, data)
732 index = None
733 try:
734 # We do not (currently) accept data for constructing Variables
--> 735 assert data is None
737 if not self.index_set().isfinite() and self._dense:
738 # Note: if the index is not finite, then we cannot
739 # iterate over it. This used to be fatal; now we
740 # just warn
741 logger.warning(
742 "Var '%s' indexed by a non-finite set, but declared "
743 "with 'dense=True'. Reverting to 'dense=False' as "
(...)
746 "'dense=False'" % (self.name,)
747 )
AssertionError:
{None: {'nodes': ['source', 'pump1', 'pump2', 'valveA', 'valveB', 'valveC'], 'pumps': ['pump1', 'pump2'], 'valves': ['valveA', 'valveB', 'valveC']}}
%% Cell type:markdown id: tags:
Frage: gibt es nur eine Lösung für Drehzahl?
Bsp. Optimierung nach Dezentraler Pumpe um modell zu prüfen
%% Cell type:code id: tags:
``` python
from pyomo.opt import SolverFactory
opt = pyo.SolverFactory('scipampl', executable=r'C:\Program Files\SCIPOptSuite 9.2.0\bin\scip.exe')
instance = modell.create_instance(TestData)
result=opt.solve(instance, tee=True)
```
%% Output
SCIP version 9.2.0 [precision: 8 byte] [memory: block] [mode: optimized] [LP solver: Soplex 7.1.2] [GitHash: 74cea9222e]
Copyright (c) 2002-2024 Zuse Institute Berlin (ZIB)
External libraries:
Soplex 7.1.2 Linear Programming Solver developed at Zuse Institute Berlin (soplex.zib.de) [GitHash: b040369c]
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)
ZIMPL 3.6.2 Zuse Institute Mathematical Programming Language developed by T. Koch (zimpl.zib.de)
AMPL/MP 690e9e7 AMPL .nl file reader library (github.com/ampl/mp)
PaPILO 2.4.0 parallel presolve for integer and linear optimization (github.com/scipopt/papilo) (built with TBB) [GitHash: 2d9fe29f]
Nauty 2.8.8 Computing Graph Automorphism Groups by Brendan D. McKay (users.cecs.anu.edu.au/~bdm/nauty)
sassy 1.1 Symmetry preprocessor by Markus Anders (github.com/markusa4/sassy)
Ipopt 3.14.16 Interior Point Optimizer developed by A. Waechter et.al. (github.com/coin-or/Ipopt)
user parameter file <scip.set> not found - using default parameters
read problem <C:\Users\STEINM~1\AppData\Local\Temp\tmptn91lkdj.pyomo.nl>
============
original problem has 5 variables (0 bin, 0 int, 0 impl, 5 cont) and 1 constraints
solve problem
=============
feasible solution found by trivial heuristic after 0.0 seconds, objective value 0.000000e+00
presolving:
(0.0s) symmetry computation started: requiring (bin +, int +, cont +), (fixed: bin -, int -, cont -)
(0.0s) symmetry computation finished: 1 generators found (max: 1500, log10 of symmetry group size: 0.0) (symcode time: 0.00)
dynamic symmetry handling statistics:
orbitopal reduction: no components
orbital reduction: no components
lexicographic reduction: 1 permutations with support sizes 4
handled 1 out of 1 symmetry components
presolving (1 rounds: 1 fast, 1 medium, 1 exhaustive):
0 deleted vars, 0 deleted constraints, 0 added constraints, 0 tightened bounds, 0 added holes, 0 changed sides, 0 changed coefficients
0 implications, 0 cliques
presolved problem has 5 variables (0 bin, 0 int, 0 impl, 5 cont) and 1 constraints
1 constraints of type <nonlinear>
Presolving Time: 0.00
transformed 3/3 original solutions to the transformed problem space
%% Cell type:code id: tags:
``` python
opt
```
%% Cell type:code id: tags:
``` python
```
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment