From 62ec7c6c1c91d24039d95632e358572172f729a1 Mon Sep 17 00:00:00 2001 From: Joris Gillis Date: Thu, 25 Jul 2013 08:21:25 +0200 Subject: [PATCH 1/5] Less verbosity for IpoptSolver --- DynOPFlow.py | 1 + 1 file changed, 1 insertion(+) diff --git a/DynOPFlow.py b/DynOPFlow.py index d267863..e77b6a5 100644 --- a/DynOPFlow.py +++ b/DynOPFlow.py @@ -708,6 +708,7 @@ def Dispatch(self, Horizon = 24, Simulation = 0, GridLoss = 'True'): solver = IpoptSolver(Cost,G) solver.setOption("parametric",True) solver.setOption("generate_hessian",True) + solver.setOption("print_level",0) solver.setOption("max_iter",2000) solver.setOption("tol",1e-4) solver.setOption("linear_solver","ma27") From 887e97b497315ef5ff7bf8261bad399fa9d8e8bb Mon Sep 17 00:00:00 2001 From: Joris Gillis Date: Thu, 25 Jul 2013 08:21:48 +0200 Subject: [PATCH 2/5] Version bumping required CasADi to 1.7.x --- DynOPFlow.py | 46 ++++++++++++++++++++++------------------------ 1 file changed, 22 insertions(+), 24 deletions(-) diff --git a/DynOPFlow.py b/DynOPFlow.py index e77b6a5..eaa3534 100644 --- a/DynOPFlow.py +++ b/DynOPFlow.py @@ -17,6 +17,8 @@ Requires the installation of the open-source Python module casADi together with the NLP solver ipopt +Required version of CasADi: v1.7.x + """ from casadi import * @@ -139,10 +141,10 @@ def setDynamics(self, RHS = [], dt = 1., nstep = 10): fimplicit = SXFunction(daeIn(x=X,p=U),daeOut(ode=RHS)) fimplicit.init() - k1 = fimplicit.eval(daeIn( x=X,p=U ))[DAE_ODE] - k2 = fimplicit.eval(daeIn( x=X+0.5*dtRK4*k1,p=U ))[DAE_ODE] - k3 = fimplicit.eval(daeIn( x=X+0.5*dtRK4*k2,p=U ))[DAE_ODE] - k4 = fimplicit.eval(daeIn( x=X+dtRK4*k3,p=U ))[DAE_ODE] + [k1] = daeOut(fimplicit.eval(daeIn( x=X,p=U )),"ode") + [k2] = daeOut(fimplicit.eval(daeIn( x=X+0.5*dtRK4*k1,p=U )),"ode") + [k3] = daeOut(fimplicit.eval(daeIn( x=X+0.5*dtRK4*k2,p=U )),"ode") + [k4] = daeOut(fimplicit.eval(daeIn( x=X+dtRK4*k3,p=U )),"ode") rk4_step = SXFunction([X,U],[X + (1./6)*dtRK4*(k1 + 2*k2 + 2*k3 + k4)]) rk4_step.init() @@ -335,15 +337,13 @@ def Flow(self, OPFSolver = 'False'): for line in range(NLine): Cost += np.real(Graph[line][2])*LineCurrents2[line] - Cost = MXFunction([V],[Cost]) - Cost.init() - - G = MXFunction([V],[g]) - G.init() - G = G.expand() + nl = MXFunction(nlpIn(V=V),nlpOut(f=Cost,g=g)) + nl.init() # set-up solver - solver = IpoptSolver(Cost,G) + solver = IpoptSolver(nl) + solver.setOption("print_level",0) + solver.setOption("expand",True) solver.setOption("parametric",False) solver.setOption("generate_hessian",True) solver.setOption("max_iter",1000) @@ -481,10 +481,10 @@ def __CostConstructor__(self, V, EP, Nstage, GridLoss): Cost += Cost_k Cost = Cost/Nstage - Cost = MXFunction([V,EP],[Cost]) - Cost.init() + CostF = MXFunction([V,EP],[Cost]) + CostF.init() - return Cost + return Cost, CostF def Dispatch(self, Horizon = 24, Simulation = 0, GridLoss = 'True'): @@ -571,12 +571,12 @@ def Dispatch(self, Horizon = 24, Simulation = 0, GridLoss = 'True'): ############################### BUILD COST AND CONSTRAINTS ############################### - Cost = self.__CostConstructor__(V, EP, Nstage, GridLoss) + Cost, CostF = self.__CostConstructor__(V, EP, Nstage, GridLoss) if (Simulation > 0): - self.CostNMPC = self.__CostConstructor__(Vstore, EP, Simulation, GridLoss) + _,self.CostNMPC = self.__CostConstructor__(Vstore, EP, Simulation, GridLoss) else: - self.CostNMPC = Cost + self.CostNMPC = CostF @@ -699,16 +699,14 @@ def Dispatch(self, Horizon = 24, Simulation = 0, GridLoss = 'True'): entry('IneqConst', expr = veccat(IneqConst)) ]) - - G = MXFunction([V,EP],[g]) - G.init() - G = G.expand() + nl = MXFunction(nlpIn(x=V,p=EP),nlpOut(f=Cost,g=g)) + nl.init() # set-up solver - solver = IpoptSolver(Cost,G) - solver.setOption("parametric",True) - solver.setOption("generate_hessian",True) + solver = IpoptSolver(nl) + solver.setOption("expand",True) solver.setOption("print_level",0) + solver.setOption("hessian_approximation","exact") solver.setOption("max_iter",2000) solver.setOption("tol",1e-4) solver.setOption("linear_solver","ma27") From 30b7b31ce97a38934ebc13c0335d9316bbf195cd Mon Sep 17 00:00:00 2001 From: Joris Gillis Date: Thu, 25 Jul 2013 09:21:35 +0200 Subject: [PATCH 3/5] Python has a boolean type: 'True' -> True, 'False' -> False --- DynOPFlow.py | 80 +++++++++++++-------------- MainOPF.py | 2 +- NMPC.py | 14 ++--- NMPCDemo.py | 4 +- NMPC_MultipleHorizon.py | 4 +- NMPC_MultipleHorizon_MonteCarlo.py | 4 +- NMPC_MultipleHorizon_MonteCarlo_II.py | 6 +- 7 files changed, 57 insertions(+), 57 deletions(-) diff --git a/DynOPFlow.py b/DynOPFlow.py index eaa3534..dae5d40 100644 --- a/DynOPFlow.py +++ b/DynOPFlow.py @@ -40,9 +40,9 @@ class Plant: - def __init__(self, Inputs = [], States = [], R = 0., Directionality = 'Mono', Load = 'False', Bus = [], label = []): # set of reasons (strings) why the Dae cannot be modified (add new x/z/u/p/output) + def __init__(self, Inputs = [], States = [], R = 0., Directionality = 'Mono', Load = False, Bus = [], label = []): # set of reasons (strings) why the Dae cannot be modified (add new x/z/u/p/output) - self._frozen = 'false' + self._frozen = False self.Bus = Bus self.label = label self.Directionality = Directionality @@ -70,7 +70,7 @@ def __init__(self, Inputs = [], States = [], R = 0., Directionality = 'Mono', Lo entry("CurrentImag")] #Structure for INPUTS of various power plants - if (Load == 'True'): + if (Load == True): InputList.append(entry("ActivePower")) InputList.append(entry("ReactivePower")) elif (Directionality == 'Bi'): @@ -110,9 +110,9 @@ def __init__(self, Inputs = [], States = [], R = 0., Directionality = 'Mono', Lo self.LB = Bound(-inf) self.UB = Bound( inf) - if (Directionality == 'Mono') and (Load == 'False'): + if (Directionality == 'Mono') and (Load == False): self.LB['Inputs','Power'] = 0. - elif (Directionality == 'Bi') and (Load == 'False'): + elif (Directionality == 'Bi') and (Load == False): self.LB['Inputs','Pcharge'] = 0. self.LB['Inputs','Pdischarge'] = 0. else: @@ -123,7 +123,7 @@ def __init__(self, Inputs = [], States = [], R = 0., Directionality = 'Mono', Lo def setDynamics(self, RHS = [], dt = 1., nstep = 10): - if (self._frozen == 'true'): + if (self._frozen == True): print "Plant already added to the grid, call ignored" return @@ -170,7 +170,7 @@ def __BuildFunc__(self, Expr, Terminal): Uprev = self.InputsPrev - if Terminal == 'False': + if Terminal == False: listFuncInput = [U, Uprev] if (X.size > 0): listFuncInput.append(X) @@ -182,9 +182,9 @@ def __BuildFunc__(self, Expr, Terminal): return Func - def setConstraints(self, Const, Terminal = 'False'): + def setConstraints(self, Const, Terminal = False): - if (self._frozen == 'true'): + if (self._frozen == True): print "Plant already added to the grid, call ignored" return @@ -193,32 +193,32 @@ def setConstraints(self, Const, Terminal = 'False'): #ConstFunc = self.__BuildFunc__(veccat(Const), Terminal) - if (Terminal == 'False'): + if (Terminal == False): self._StageConst = self.__BuildFunc__(veccat(Const), Terminal) - elif (Terminal == 'True'): + elif (Terminal == True): self._TerminalConst = self.__BuildFunc__(veccat(Const), Terminal) - def setCost(self, Cost, Terminal = 'False'): - if (self._frozen == 'true'): + def setCost(self, Cost, Terminal = False): + if (self._frozen == True): print "Plant already added to the grid, call ignored" return #CostFunc = self.__BuildFunc__(Cost, Terminal) - if (Terminal == 'False'): + if (Terminal == False): self._StageCost = self.__BuildFunc__(Cost, Terminal) - elif (Terminal == 'True'): + elif (Terminal == True): self._TerminalCost = self.__BuildFunc__(Cost, Terminal) def addPlant(self, Net): - if (self._frozen == 'true'): + if (self._frozen == True): print "Plant already added to the grid, call ignored" return Net.PlantList.append(self) - self._frozen = 'true' + self._frozen = True class PowerGrid: """ @@ -239,7 +239,7 @@ def __init__(self, NBus = 0,Graph = []): #CONSTRUCT POWER FLOW - def Flow(self, OPFSolver = 'False'): + def Flow(self, OPFSolver = False): NBus = self.NBus NLine = np.size( self.Graph ,axis = 0) @@ -302,10 +302,10 @@ def Flow(self, OPFSolver = 'False'): self.BusCurrentsImagFunc.init() self.LineCurrents2Func.init() self.BusVoltages2Func.init() - self.OPF = 'True' + self.OPF = True #CONSTRUCT OPF SOLVER IS ASKED IN THE FLOW FUNCTION OPTIONS - if (OPFSolver == 'True'): + if (OPFSolver == True): print "Construct OPF Solver" Power = struct_ssym([ @@ -370,12 +370,12 @@ def OPFSolve(self, Grid = []): if not(hasattr(self,'OPF')): #Check that .Flow() has been called - print "You must call .Flow(OPFSolver = 'True') to setup OPF before calling .OPFSolve()" + print "You must call .Flow(OPFSolver = True) to setup OPF before calling .OPFSolve()" return [] - if (self.OPF == 'True'): + if (self.OPF == True): #Check that a solver exists - print "You must call .Flow(OPFSolver = 'True') to setup OPF before calling .OPFSolve()" + print "You must call .Flow(OPFSolver = True) to setup OPF before calling .OPFSolve()" return [] @@ -453,7 +453,7 @@ def __CostConstructor__(self, V, EP, Nstage, GridLoss): Cost = 0 - if (GridLoss == 'True'): + if (GridLoss == True): NLine = len( self.Graph ) for k in range(Nstage): # Grid losses @@ -487,11 +487,11 @@ def __CostConstructor__(self, V, EP, Nstage, GridLoss): return Cost, CostF - def Dispatch(self, Horizon = 24, Simulation = 0, GridLoss = 'True'): + def Dispatch(self, Horizon = 24, Simulation = 0, GridLoss = True): """ Constructs the power dispatch problem, default Horizon length (if argument Horizon is not provided) is 24 time units """ - if (self.OPF == 'False'): + if (self.OPF == False): Power.Flow(self) print "Construct Dynamic OPF" @@ -634,7 +634,7 @@ def Dispatch(self, Horizon = 24, Simulation = 0, GridLoss = 'True'): # Plant participating current squared, i.e. |i|**2 PlantCurrent2 = PlantCurrentReal*PlantCurrentReal + PlantCurrentImag*PlantCurrentImag - if (plant._Load == 'True'): + if (plant._Load == True): #Load fixing: [Active, Reactive] = Consumed Active / Reactive power EquConst.append(ParticipatingActivePower - V['Inputs',k,plant.label,'ActivePower']) EquConst.append(ParticipatingReactivePower - V['Inputs',k,plant.label,'ReactivePower']) @@ -799,7 +799,7 @@ def Profiles(self, N = 0): #ASSIGN PROFILES & SOLVE - def DYNSolve(self, x0 = [], u0 = 0., init = [], time = 0, Periodic = 'False', EW = 1.): + def DYNSolve(self, x0 = [], u0 = 0., init = [], time = 0, Periodic = False, EW = 1.): lbV = self.VOptDispatch(-inf) ubV = self.VOptDispatch( inf) @@ -851,7 +851,7 @@ def DYNSolve(self, x0 = [], u0 = 0., init = [], time = 0, Periodic = 'False', EW lbV.cat = lbV.cat - (ubV.cat == inf)*ubV.cat #Set inf in x0 to free the initial conditions ###### PERIODIC CONSTRAINTS (IF REQUIRED) ####### - #if (Periodic == 'False'): + #if (Periodic == False): # lbg['Periodic'] = -inf # ubg['Periodic'] = inf @@ -931,7 +931,7 @@ def NMPCSimulation(self, x0 = [], u0 = [], init = [], Simulation = 0): #Extract results - def ExtractInfo(self, v_opt, BusPower = 'False', PlantPower = 'False', TotalPower = 'False' ): + def ExtractInfo(self, v_opt, BusPower = False, PlantPower = False, TotalPower = False ): self.SolutionInfo = {} @@ -961,15 +961,15 @@ def ExtractInfo(self, v_opt, BusPower = 'False', PlantPower = 'False', TotalPowe self.SolutionInfo['LineCurrentsModule'] = np.concatenate(LineCurrents_opt,axis=1) #### Total Powers - if (TotalPower == 'True'): + if (TotalPower == True): TotalPower = {} TotalPower['Load'] = 0 - for plant in [plant for plant in self.PlantList if (plant._Load == 'True')]: + for plant in [plant for plant in self.PlantList if (plant._Load == True)]: TotalPower['Load'] -= np.array(v_opt['Inputs',:,plant.label,'ActivePower']) TotalPower['Injected'] = 0. - for plant in [plant for plant in self.PlantList if not(plant._Load == 'True')]: + for plant in [plant for plant in self.PlantList if not(plant._Load == True)]: if (plant.Directionality == 'Mono'): TotalPower['Injected'] += np.array(v_opt['Inputs',:,plant.label,'Power']) else: @@ -979,7 +979,7 @@ def ExtractInfo(self, v_opt, BusPower = 'False', PlantPower = 'False', TotalPowe #Construct implicit values in the network - if (BusPower == 'True'): + if (BusPower == True): BusActivePower = [] BusReactivePower = [] @@ -1007,7 +1007,7 @@ def ExtractInfo(self, v_opt, BusPower = 'False', PlantPower = 'False', TotalPowe self.SolutionInfo['BusCurrentModule'] = np.concatenate( BusCurrentModule, axis=0) self.SolutionInfo['BusCurrentAngle'] = np.concatenate( BusCurrentAngle, axis=0) - if (PlantPower == 'True'): + if (PlantPower == True): PlantActivePowerDictionary = {} PlantReactivePowerDictionary = {} CosPhiDictionary = {} @@ -1049,7 +1049,7 @@ def ExtractInfo(self, v_opt, BusPower = 'False', PlantPower = 'False', TotalPowe #### RESULT PLOTTING ##### - def DYNSolvePlot(self, v_opt, NMPC = 'False', dt = 1, Path = []): + def DYNSolvePlot(self, v_opt, NMPC = False, dt = 1, Path = []): SavedFigs = [] @@ -1081,7 +1081,7 @@ def SaveFig(Path,Name): time[key] = np.array([k*dt for k in range(len(v_opt[key]))]).T # construct a list of the plants (excluding the loads) - PlantList = [plant for plant in self.PlantList if not(plant._Load == 'True')] + PlantList = [plant for plant in self.PlantList if not(plant._Load == True)] SizeSubplt = np.ceil(sqrt(len(PlantList))) SizeSubpltAll = np.ceil(sqrt(len(self.PlantList))) @@ -1166,7 +1166,7 @@ def SaveFig(Path,Name): plt.hold('on') for plant in self.PlantList: Power = [] - if plant._Load == 'True': + if plant._Load == True: Power = -np.array(v_opt['Inputs',:,plant.label,'ActivePower']) elif plant.Directionality == 'Bi': Power = np.array(v_opt['Inputs',:,plant.label,'Pdischarge']) - np.array(v_opt['Inputs',:,plant.label,'Pcharge']) @@ -1235,7 +1235,7 @@ def SaveFig(Path,Name): plt.show() return SavedFigs - def DYNSolvePlotCompare(self, v_opt, NMPC = 'False', label = '', marker = '', linewidth = 1): + def DYNSolvePlotCompare(self, v_opt, NMPC = False, label = '', marker = '', linewidth = 1): Nstage = self.TimeSetup['Horizon'] @@ -1310,4 +1310,4 @@ def DYNSolvePlotCompare(self, v_opt, NMPC = 'False', label = '', marker = '', li plt.title('Wind Speed m/s') - \ No newline at end of file + diff --git a/MainOPF.py b/MainOPF.py index fa82662..b4d6d4f 100644 --- a/MainOPF.py +++ b/MainOPF.py @@ -32,7 +32,7 @@ -Net.Flow(OPFSolver = 'True') +Net.Flow(OPFSolver = True) diff --git a/NMPC.py b/NMPC.py index d0ad4a3..2078111 100644 --- a/NMPC.py +++ b/NMPC.py @@ -69,9 +69,9 @@ Hydro.setCost ( Cost ) #TConst = (15-h) -#Hydro.setConstraints(TConst, Terminal = 'True') +#Hydro.setConstraints(TConst, Terminal = True) #TCost = -100*h -#Hydro.setCost(TCost, Terminal = 'True') +#Hydro.setCost(TCost, Terminal = True) Hydro.addPlant(Net) @@ -157,7 +157,7 @@ Thermal.UB['Inputs','Power'] = 1000 ##### Load ###### -Load = Plant(Load = 'True', Bus = 0, label = 'Load') +Load = Plant(Load = True, Bus = 0, label = 'Load') Load.addPlant(Net) Load.LB['Inputs', 'ActivePower'] = -1000 Load.LB['Inputs','ReactivePower'] = -750 @@ -224,16 +224,16 @@ Sol, stats = Net.DYNSolve(x0 = x0, u0 = u0, time = 0, init = init) init = Sol -Net.ExtractInfo(Sol, PlantPower = 'True', BusPower = 'True', TotalPower = 'True') -Net.DYNSolvePlot(Sol, dt = 1) +Net.ExtractInfo(Sol, PlantPower = True, BusPower = True, TotalPower = True) +#Net.DYNSolvePlot(Sol, dt = 1) Traj, NMPC_Info = Net.NMPCSimulation(x0 = x0, u0 = u0, init = init, Simulation = Nsimulation) #Exctract info -Net.ExtractInfo(Traj, PlantPower = 'True', BusPower = 'True', TotalPower = 'True') -Net.DYNSolvePlot(Traj, dt = 1) +Net.ExtractInfo(Traj, PlantPower = True, BusPower = True, TotalPower = True) +#Net.DYNSolvePlot(Traj, dt = 1) diff --git a/NMPCDemo.py b/NMPCDemo.py index 9a0f4a3..0487a3d 100644 --- a/NMPCDemo.py +++ b/NMPCDemo.py @@ -132,7 +132,7 @@ Thermal.UB['Inputs','Power'] = 1000 ##### Load ###### -Load = Plant(Load = 'True', Bus = 0, label = 'Load') +Load = Plant(Load = True, Bus = 0, label = 'Load') Load.addPlant(Net) Load.LB['Inputs', 'ActivePower'] = -1000 Load.LB['Inputs','ReactivePower'] = -750 @@ -184,7 +184,7 @@ Traj, NMPC_Info = Net.NMPCSimulation(x0 = x0, u0 = u0, init = init, Simulation = Nsimulation) #Plotting -Net.ExtractInfo(Traj, PlantPower = 'True', BusPower = 'True', TotalPower = 'True') +Net.ExtractInfo(Traj, PlantPower = True, BusPower = True, TotalPower = True) Path = '/Users/sebastien/Desktop/Research/PowerFlowCodes/Paper/Figures' SavedFigs = Net.DYNSolvePlot(Traj, dt = 1/24., Path = Path) diff --git a/NMPC_MultipleHorizon.py b/NMPC_MultipleHorizon.py index 7c048bd..92b2626 100644 --- a/NMPC_MultipleHorizon.py +++ b/NMPC_MultipleHorizon.py @@ -146,7 +146,7 @@ Thermal.UB['Inputs','Power'] = 1000 ##### Load ###### -Load = Plant(Load = 'True', Bus = 0, label = 'Load') +Load = Plant(Load = True, Bus = 0, label = 'Load') Load.addPlant(Net) Load.LB['Inputs', 'ActivePower'] = -1000 Load.LB['Inputs','ReactivePower'] = -750 @@ -222,7 +222,7 @@ plt.plot(Horizons,Cost,linestyle = 'none',marker = 'o') plt.show() #Exctract info -#Net.ExtractInfo(Net.Vstore, PlantPower = 'True', BusPower = 'True', TotalPower = 'True') +#Net.ExtractInfo(Net.Vstore, PlantPower = True, BusPower = True, TotalPower = True) #Net.DYNSolvePlot(Net.Vstore, dt = 1) diff --git a/NMPC_MultipleHorizon_MonteCarlo.py b/NMPC_MultipleHorizon_MonteCarlo.py index a00a939..2b52fb4 100644 --- a/NMPC_MultipleHorizon_MonteCarlo.py +++ b/NMPC_MultipleHorizon_MonteCarlo.py @@ -147,7 +147,7 @@ Thermal.UB['Inputs','Power'] = 1000 ##### Load ###### -Load = Plant(Load = 'True', Bus = 0, label = 'Load') +Load = Plant(Load = True, Bus = 0, label = 'Load') Load.addPlant(Net) Load.LB['Inputs', 'ActivePower'] = -1000 Load.LB['Inputs','ReactivePower'] = -750 @@ -236,7 +236,7 @@ plt.plot(Horizons,Cost,linestyle = 'none',marker = 'o') plt.show() #Exctract info -#Net.ExtractInfo(Net.Vstore, PlantPower = 'True', BusPower = 'True', TotalPower = 'True') +#Net.ExtractInfo(Net.Vstore, PlantPower = True, BusPower = True, TotalPower = True) #Net.DYNSolvePlot(Net.Vstore, dt = 1) diff --git a/NMPC_MultipleHorizon_MonteCarlo_II.py b/NMPC_MultipleHorizon_MonteCarlo_II.py index 18280d3..3ee9c96 100644 --- a/NMPC_MultipleHorizon_MonteCarlo_II.py +++ b/NMPC_MultipleHorizon_MonteCarlo_II.py @@ -164,7 +164,7 @@ Thermal.UB['Inputs','Power'] = 1000 ##### Load ###### -Load = Plant(Load = 'True', Bus = 0, label = 'Load') +Load = Plant(Load = True, Bus = 0, label = 'Load') Load.addPlant(Net) Load.LB['Inputs', 'ActivePower'] = -1000 Load.LB['Inputs','ReactivePower'] = -750 @@ -238,7 +238,7 @@ Sol, stats = Net.DYNSolve(x0 = x0, u0 = u0, time = 0, init = init) init = Sol - Net.ExtractInfo(Sol, PlantPower = 'True', BusPower = 'True', TotalPower = 'True') + Net.ExtractInfo(Sol, PlantPower = True, BusPower = True, TotalPower = True) Net.DYNSolvePlot(Sol, dt = 1) assert(0==1) ##### NMPC Simulations ##### @@ -250,7 +250,7 @@ else: Results[-1]['Cost'] += NMPC_Info['Cost']/float(HorizonSetup['Nrealization']) - #Net.ExtractInfo(Traj, PlantPower = 'True', BusPower = 'True', TotalPower = 'True') + #Net.ExtractInfo(Traj, PlantPower = True, BusPower = True, TotalPower = True) #Net.DYNSolvePlot(Traj, dt = 1) Horizon += HorizonSetup['Step'] From d83774672429cb2e89a7314065042eaca3e6e793 Mon Sep 17 00:00:00 2001 From: Joris Gillis Date: Thu, 25 Jul 2013 09:46:21 +0200 Subject: [PATCH 4/5] python-style correction: under-under names should not be contructed --- DynOPFlow.py | 20 ++++++++++---------- 1 file changed, 10 insertions(+), 10 deletions(-) diff --git a/DynOPFlow.py b/DynOPFlow.py index dae5d40..5bf8d3a 100644 --- a/DynOPFlow.py +++ b/DynOPFlow.py @@ -163,7 +163,7 @@ def setDynamics(self, RHS = [], dt = 1., nstep = 10): self._Shoot = Shoot - def __BuildFunc__(self, Expr, Terminal): + def _BuildFunc(self, Expr, Terminal): X = self.States U = self.Inputs @@ -191,12 +191,12 @@ def setConstraints(self, Const, Terminal = False): if not(isinstance(Const,list)): Const = [Const] - #ConstFunc = self.__BuildFunc__(veccat(Const), Terminal) + #ConstFunc = self._BuildFunc(veccat(Const), Terminal) if (Terminal == False): - self._StageConst = self.__BuildFunc__(veccat(Const), Terminal) + self._StageConst = self._BuildFunc(veccat(Const), Terminal) elif (Terminal == True): - self._TerminalConst = self.__BuildFunc__(veccat(Const), Terminal) + self._TerminalConst = self._BuildFunc(veccat(Const), Terminal) def setCost(self, Cost, Terminal = False): @@ -204,12 +204,12 @@ def setCost(self, Cost, Terminal = False): print "Plant already added to the grid, call ignored" return - #CostFunc = self.__BuildFunc__(Cost, Terminal) + #CostFunc = self._BuildFunc(Cost, Terminal) if (Terminal == False): - self._StageCost = self.__BuildFunc__(Cost, Terminal) + self._StageCost = self._BuildFunc(Cost, Terminal) elif (Terminal == True): - self._TerminalCost = self.__BuildFunc__(Cost, Terminal) + self._TerminalCost = self._BuildFunc(Cost, Terminal) def addPlant(self, Net): @@ -446,7 +446,7 @@ def OPFSolve(self, Grid = []): ########### POWER DISPACTH PROBLEM ########## - def __CostConstructor__(self, V, EP, Nstage, GridLoss): + def _CostConstructor(self, V, EP, Nstage, GridLoss): """ Constructor for the Cost function, handy to build the cost for different V:s """ @@ -571,10 +571,10 @@ def Dispatch(self, Horizon = 24, Simulation = 0, GridLoss = True): ############################### BUILD COST AND CONSTRAINTS ############################### - Cost, CostF = self.__CostConstructor__(V, EP, Nstage, GridLoss) + Cost, CostF = self._CostConstructor(V, EP, Nstage, GridLoss) if (Simulation > 0): - _,self.CostNMPC = self.__CostConstructor__(Vstore, EP, Simulation, GridLoss) + _,self.CostNMPC = self._CostConstructor(Vstore, EP, Simulation, GridLoss) else: self.CostNMPC = CostF From d324705b96b54821b9322231ebc2edfc5f33dbad Mon Sep 17 00:00:00 2001 From: Joris Gillis Date: Thu, 25 Jul 2013 16:57:31 +0200 Subject: [PATCH 5/5] Added automatic scaling of optimiation problem --- DynOPFlow.py | 109 +++++++++++++++++++++++++++++++++++++++++++++------ NMPCDemo.py | 2 +- 2 files changed, 97 insertions(+), 14 deletions(-) diff --git a/DynOPFlow.py b/DynOPFlow.py index b30f820..4cb65c7 100644 --- a/DynOPFlow.py +++ b/DynOPFlow.py @@ -36,8 +36,23 @@ #rho_air = 1.2 #rho_water = 1e3 #gravity = 9.81 - - + +# Helper function to obtain automatic scaling +class BadScalingException(Exception): + pass + +scale_min = 1e-5 + +def scaleFromBounds(lb,ub): + lb = filter(lambda x: x!=-inf,lb) + ub = filter(lambda x: x!=inf,ub) + lb = min(lb) + ub = max(ub) + scale = max(abs(lb),abs(ub)) + if scale < scale_min: + raise BadScalingException() + else: + return scale class Plant: def __init__(self, Inputs = [], States = [], R = 0., Directionality = 'Mono', Load = False, Bus = [], label = []): # set of reasons (strings) why the Dae cannot be modified (add new x/z/u/p/output) @@ -696,7 +711,7 @@ def Dispatch(self, Horizon = 24, Simulation = 0, GridLoss = True): # set-up solver solver = IpoptSolver(nl) solver.setOption("expand",True) - solver.setOption("print_level",0) + #solver.setOption("print_level",0) solver.setOption("hessian_approximation","exact") solver.setOption("max_iter",2000) solver.setOption("tol",1e-4) @@ -715,6 +730,11 @@ def Dispatch(self, Horizon = 24, Simulation = 0, GridLoss = True): self.OptDispatch = solver self.gOptDispatch = g + self.scaling_V = self.VOptDispatch(1) + self.scaling_g = self.gOptDispatch(1) + self.scaling_f = 1 + self.OptDispatch_nl = nl + self.Properties = BusProperties print self.Properties ############## SOLVER CONSTRUCTED ############## @@ -790,7 +810,7 @@ def Profiles(self, N = 0): #ASSIGN PROFILES & SOLVE - def DYNSolve(self, x0 = [], u0 = 0., init = [], time = 0, Periodic = False, EW = 1.): + def DYNSolve(self, x0 = [], u0 = 0., init = [], time = 0, Periodic = False, EW = 1., resetScale=False): lbV = self.VOptDispatch(-inf) ubV = self.VOptDispatch( inf) @@ -844,7 +864,7 @@ def DYNSolve(self, x0 = [], u0 = 0., init = [], time = 0, Periodic = False, EW = lbV['States',0] = x0 ubV['States',0] = x0 - lbV.cat = lbV.cat - (ubV.cat == inf)*ubV.cat #Set inf in x0 to free the initial conditions + #lbV.cat = lbV.cat - (ubV.cat == inf)*ubV.cat #Set inf in x0 to free the initial conditions ###### PERIODIC CONSTRAINTS (IF REQUIRED) ####### #if (Periodic == False): @@ -854,21 +874,82 @@ def DYNSolve(self, x0 = [], u0 = 0., init = [], time = 0, Periodic = False, EW = EP = self.EP() EP['u0'] = u0 EP['EW'] = EW + + + # Logic to calculate scaling + if resetScale: + + # Scaling for decision variables + self.scaling_V = self.VOptDispatch(1) + for k in ["Inputs","States"]: + for c in self.VOptDispatch.getStruct(k).canonicalIndices(): + i = (k,slice(None))+c + try: + scale = self.scaling_V.__setitem__(i,scaleFromBounds(lbV.__getitem__(i),ubV.__getitem__(i))) + except BadScalingException: + print "Warning. Variable %s is badly scaled by bounds." % str(i) + for bus in range(NBus): + scale = self.PowerFlowBounds['Vmax'][bus] + if scale < scale_min or scale==inf: + print "Warning. Bus VMax %d is badly scaled." % str(bus) + else: + self.scaling_V["BusVoltages",:,"Imag",bus] = scale + self.scaling_V["BusVoltages",:,"Real",bus] = scale + + # Scaling for nonlinear bounds + self.scaling_g = self.gOptDispatch(1) + for k in ["CurrentBalance","BusVoltages2","LineCurrents2"]: + for i in range(lbg[k,0].size()): + try: + self.scaling_g[k,:,i] = scaleFromBounds(lbg[k,:,i],ubg[k,:,i]) + except: + print "Warning. Constraint (%s,:,%d) is badly scaled by bounds." % (k,i) + + # Scaling for objective function + nl = self.OptDispatch_nl + X, P = nlpIn(nl.symbolicInput(),"x","p") + + nl.setInput(EP, "p") + nl.setInput(lbV, "x") + nl.evaluate() + f_low = nl.getOutput()/10 + nl.setInput(ubV, "x") + nl.evaluate() + f_high = nl.getOutput()/10 + nl.setInput(init, "x") + nl.evaluate() + f_init = nl.getOutput() + + self.scaling_f = max([abs(f_low),abs(f_high),abs(f_init)]) + + # Construct an alternative solver + F,G = nlpOut(nl.call(nlpIn(x=X*self.scaling_V.cat,p=P)),"f","g") + + scaled_nl = MXFunction(nlpIn(x=X,p=P),nlpOut(f=F/self.scaling_f,g=G/self.scaling_g.cat)) + + options = self.OptDispatch.dictionary() + self.OptDispatch = IpoptSolver(scaled_nl) + self.OptDispatch.setOption(options) + self.OptDispatch.init() + + - self.OptDispatch.setInput(lbV, "lbx") - self.OptDispatch.setInput(ubV, "ubx") - self.OptDispatch.setInput(init, "x0" ) - self.OptDispatch.setInput(lbg, "lbg") - self.OptDispatch.setInput(ubg, "ubg") + self.OptDispatch.setInput(lbV/self.scaling_V.cat, "lbx") + self.OptDispatch.setInput(ubV/self.scaling_V.cat, "ubx") + self.OptDispatch.setInput(init/self.scaling_V.cat, "x0" ) + self.OptDispatch.setInput(lbg/self.scaling_g.cat, "lbg") + self.OptDispatch.setInput(ubg/self.scaling_g.cat, "ubg") self.OptDispatch.setInput(EP, "p") + + self.OptDispatch.solve() self.lbV = lbV self.ubV = ubV - v_opt = self.VOptDispatch(self.OptDispatch.output("x")) + v_opt = self.VOptDispatch(self.OptDispatch.output("x")*self.scaling_V.cat) success = int(self.OptDispatch.getStat('return_status') == 'Solve_Succeeded') @@ -899,13 +980,15 @@ def Simulate(self, Sol, x0, u0): - def NMPCSimulation(self, x0 = [], u0 = [], init = [], Simulation = 0): + def NMPCSimulation(self, x0 = [], u0 = [], init = [], Simulation = 0,resetScale=False): ##### NMPC Loop ##### NMPC = {'time': 0, 'success' : [], 'Traj' : []} Vstore = self.Vstore while (NMPC['time'] < Simulation): - Sol, stats = self.DYNSolve(x0 = x0, u0 = u0, time = NMPC['time'], init = init) + Sol, stats = self.DYNSolve(x0 = x0, u0 = u0, time = NMPC['time'], init = init, resetScale=resetScale) + + resetScale = False NMPC['success'].append(stats) NMPC['Traj'].append(Sol) diff --git a/NMPCDemo.py b/NMPCDemo.py index 652be1f..a5fea5f 100644 --- a/NMPCDemo.py +++ b/NMPCDemo.py @@ -213,7 +213,7 @@ #Net.DYNSolvePlot(Sol, dt = 1) # #assert(0==1) -Traj, NMPC_Info = Net.NMPCSimulation(x0 = x0, u0 = u0, init = init, Simulation = Nsimulation) +Traj, NMPC_Info = Net.NMPCSimulation(x0 = x0, u0 = u0, init = init, Simulation = Nsimulation, resetScale=True) #Plotting Net.ExtractInfo(Traj, PlantPower = True, BusPower = True, TotalPower = True)