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)