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

rebuild continuity constraint after deletion

parent ad21c566
No related branches found
No related tags found
No related merge requests found
Show whitespace changes
Inline Side-by-side
Showing
with 12 additions and 54 deletions
optimierung_pymoo.ipynb
+ 12
− 54
View file @ b8fb9bf2
%% 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,i} \geq \sum_{strang} Q_v + \sum_{strang} Q_p \\ Q_{p,i} \geq \sum_{strang} Q_v + \sum_{strang} Q_p \\
Q_p , n\epsilon [n_{min},n_{max}] \\ Q_p , 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
``` ```
%% Output %% Output
Defaulting to user installation because normal site-packages is not writeable Defaulting to user installation because normal site-packages is not writeable
Requirement already satisfied: pyomo in c:\users\steinmann\appdata\roaming\python\python312\site-packages (6.8.2) Requirement already satisfied: pyomo in c:\users\steinmann\appdata\roaming\python\python312\site-packages (6.8.2)
Requirement already satisfied: ply in c:\users\steinmann\appdata\roaming\python\python312\site-packages (from pyomo) (3.11) Requirement already satisfied: ply in c:\users\steinmann\appdata\roaming\python\python312\site-packages (from pyomo) (3.11)
[notice] A new release of pip is available: 24.2 -> 25.0 [notice] A new release of pip is available: 24.2 -> 25.0
[notice] To update, run: C:\Program Files\Python312\python.exe -m pip install --upgrade pip [notice] To update, run: C:\Program Files\Python312\python.exe -m pip install --upgrade pip
%% 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)
``` ```
%% 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 löst sich auf wenn die nachfolgenden Durchflüsse klar sind mit $-l \sum_{i \in \mathrm{S}} Q^2_i = \alpha_1 Q^2 + \alpha_2 Q n + \alpha_3 n^2$ 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) modell.Q_valve=pyo.Param(modell.valves)
#Optimierungsvariable #Optimierungsvariable
modell.n = pyo.Var(modell.pumps,bounds=(750/3600,1)) modell.n = pyo.Var(modell.pumps,bounds=(750/3600,1))
modell.Q = pyo.Var(modell.nodes) modell.Q = pyo.Var(modell.nodes)
#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[i]**3,(modell.Q[i]**2)*modell.n[i],modell.Q[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): def PumpFlow(modell,pump):
return np.dot(np.array([modell.Q[pump]**2,modell.n[pump]*modell.Q[pump],modell.n[pump]**2]),LR_H.coef_) pump=np.dot(np.array([modell.Q[pump]**2,modell.n[pump]*modell.Q[pump],modell.n[pump]**2]),LR_H.coef_)
def Flow_Equation(modell,p): return pump>=sum(modell.Q_valve[node] for node in graph.successors(pump) if node in modell.valves)+sum(modell.Q[n] for n in graph.successors(node) if node in modell.pumps)
return PumpFlow(modell=modell,pump=p) - pyo.summation(modell.Q,index=graph.successors(node)) #modell.Flow_Objective = pyo.Objective(rule=PumpFlow,sense=pyo.as_boolean)
modell.Flow_Objective = pyo.Objective(modell.pumps,rule=Flow_Equation,sense=pyo.minimize)
#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[i] for i in graph.predecessors(node))==sum(modell.Q[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
print([key for key in graph.nodes.keys() if 'valve' in key])
TestData={ TestData={
None:{ None:{
'nodes':[key for key in graph.nodes.keys()], 'nodes':[key for key in graph.nodes.keys()],
'pumps':[key for key in graph.nodes.keys() if 'pump' 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], 'valves':[key for key in graph.nodes.keys() if 'valve' in key],
'Q_valve':{'valveA':4.,'valveB':4.,'valveC':4.} 'Q_valve':{'valveA':4.,'valveB':4.,'valveC':4.},
} }
} }
print(TestData) print(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.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[i],modell.n[i],LR_P) for i in modell.pumps),sense=min)
``` ```
%% Output %% Output
['valveA', 'valveB', 'valveC']
{None: {'nodes': ['source', 'pump1', 'pump2', 'valveA', 'valveB', 'valveC'], 'pumps': ['pump1', 'pump2'], 'valves': ['valveA', 'valveB', 'valveC'], 'Q_valve': {'valveA': 4.0, 'valveB': 4.0, 'valveC': 4.0}}} {None: {'nodes': ['source', 'pump1', 'pump2', 'valveA', 'valveB', 'valveC'], 'pumps': ['pump1', 'pump2'], 'valves': ['valveA', 'valveB', 'valveC'], 'Q_valve': {'valveA': 4.0, 'valveB': 4.0, 'valveC': 4.0}}}
%% 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)
print(instance) instance.Continuity_constaint=pyo.Constraint(instance.nodes, rule=continuityRule)
instance.continuity_check=pyo.Constraint(instance.nodes,rule=continuityRule)
#result=opt.solve(instance, tee=True) #result=opt.solve(instance, tee=True)
``` ```
%% Output
unknown
invalid parameter <C:\Users\STEINM~1\AppData\Local\Temp\tmp1r989oiw.pyomo.nl>
invalid parameter <-AMPL>
syntax: C:\Program Files\SCIPOptSuite 7.0.2\bin\scip.exe [-l <logfile>] [-q] [-s <settings>] [-r <randseed>] [-f <problem>] [-b <batchfile>] [-c "command"]
-v, --version : print version and build options
-l <logfile> : copy output into log file
-q : suppress screen messages
-s <settings> : load parameter settings (.set) file
-f <problem> : load and solve problem file
-o <primref> <dualref> : pass primal and dual objective reference values for validation at the end of the solve
-b <batchfile>: load and execute dialog command batch file (can be used multiple times)
-r <randseed> : nonnegative integer to be used as random seed. Has priority over random seed specified through parameter settings (.set) file
-c "command" : execute single line of dialog commands (can be used multiple times)
---------------------------------------------------------------------------
FileNotFoundError Traceback (most recent call last)
Cell In[46], line 8
5 print(instance)
6 instance.continuity_check=pyo.Constraint(instance.nodes,rule=continuityRule)
----> 8 result=opt.solve(instance, tee=True)
File ~\AppData\Roaming\Python\Python312\site-packages\pyomo\opt\base\solvers.py:636, in OptSolver.solve(self, *args, **kwds)
630 if self._report_timing:
631 print(
632 " %6.2f seconds required for solver"
633 % (solve_completion_time - presolve_completion_time)
634 )
--> 636 result = self._postsolve()
637 result._smap_id = self._smap_id
638 result._smap = None
File ~\AppData\Roaming\Python\Python312\site-packages\pyomo\solvers\plugins\solvers\SCIPAMPL.py:222, in SCIPAMPL._postsolve(self)
216 version = self._known_versions[executable]
218 if version < (8, 0, 0, 0):
219 # it may be possible to get results from older version but this was
220 # not tested, so the old way of doing things is here preserved
--> 222 results = super(SCIPAMPL, self)._postsolve()
224 else:
225 # repeat code from super(SCIPAMPL, self)._postsolve()
226 # in order to access the log file and get the results from there
228 if self._log_file is not None:
File ~\AppData\Roaming\Python\Python312\site-packages\pyomo\opt\solver\shellcmd.py:290, in SystemCallSolver._postsolve(self)
287 results = None
289 if self._results_format is not None:
--> 290 results = self.process_output(self._rc)
291 #
292 # If keepfiles is true, then we pop the
293 # TempfileManager context while telling it to
294 # _not_ remove the files.
295 #
296 if not self._keepfiles:
297 # in some cases, the solution filename is
298 # not generated via the temp-file mechanism,
299 # instead being automatically derived from
300 # the input lp/nl filename. so, we may have
301 # to clean it up manually.
File ~\AppData\Roaming\Python\Python312\site-packages\pyomo\opt\solver\shellcmd.py:395, in SystemCallSolver.process_output(self, rc)
388 results = self._results_reader(
389 self._results_file,
390 res=results,
391 soln=results.solution(0),
392 suffixes=self._suffixes,
393 )
394 else:
--> 395 results = self._results_reader(
396 self._results_file, res=results, suffixes=self._suffixes
397 )
398 results_reader_completion_time = time.time()
399 if self._report_timing is True:
File ~\AppData\Roaming\Python\Python312\site-packages\pyomo\opt\plugins\sol.py:40, in ResultsReader_sol.__call__(self, filename, res, soln, suffixes)
36 """
37 Parse a ``*.sol`` file
38 """
39 try:
---> 40 with open(filename, "r") as f:
41 return self._load(f, res, soln, suffixes)
42 except ValueError as e:
FileNotFoundError: [Errno 2] No such file or directory: 'C:\\Users\\STEINM~1\\AppData\\Local\\Temp\\tmp1r989oiw.pyomo.sol'
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment