#Plot_q_funcs.py #Solves the optimal stopping problem as outlined in Figure 1, produces Figure 2 based on q_funcs.npy data import numpy as np import matplotlib.pyplot as plt fig = plt.figure(figsize=(8,12)) plt.rc('font', family='serif') plt.rc('font', serif='Times New Roman') plt.rc('font', size=14) q_funcs=np.load('q_funcs.npy') x=q_funcs[:,0] F_proc_A=q_funcs[:,1] F_proc_B=q_funcs[:,2] E_proc_A=q_funcs[:,3] E_proc_B=q_funcs[:,4] # Valid filling interval 1 [20161.0,28123.0] index [0,7962] # Valid emptying interval 1 [28124.0.0,32570.0] index [7963,12409] # Valid filling interval 2 [32571.0, 38357.0] index [12410,18196] # Valid emptying interval 2 [38358.0,54087.0] index [18197,33926] #Construct l_A and l_B for valid filling interval 1 and valid emptying interval 1 #This proves unsuccessful #Initialise arrays a_min_id=0 a_max_id=7962 b_min_id=7963 b_max_id=12409 array_size = a_max_id-a_min_id l_A11=np.zeros((array_size,2)) l_B11=np.zeros((array_size,2)) #Run through possible charging points to construct l_A for j in range(array_size): i=a_min_id+j l_A11[j,0]=x[i] val_u=F_proc_A[i] #Does val_u lie in the range of possible values in the discharging region if val_u<=np.nanmax(E_proc_A[b_min_id:b_max_id]) and val_u>=np.nanmin(E_proc_A[b_min_id:b_max_id]): #find the value of l_A(u) k=np.where(E_proc_A[b_min_id:b_max_id]<=val_u)[0][0] val_delta=val_u-E_proc_A[b_min_id+k-1] int_delta= E_proc_A[b_min_id+k]-E_proc_A[b_min_id+k-1] ratio=val_delta/int_delta l_A11[j,1]=x[b_min_id+k-1]+ratio*(x[b_min_id+k]-x[b_min_id+k-1]) else: #No match, ignore l_A11[j,1]=np.nan #Run through possible charging points to construct l_B for j in range(array_size): i=a_min_id+j l_B11[j,0]=x[i] val_u=F_proc_B[i] #Does val_u lie in the range of possible values in the discharging region if val_u<=np.nanmax(E_proc_B[b_min_id:b_max_id]) and val_u>=np.nanmin(E_proc_B[b_min_id:b_max_id]): #find the value of l_B(u) k=np.where(E_proc_B[b_min_id:b_max_id]<=val_u)[0][0] val_delta=val_u-E_proc_B[b_min_id+k] int_delta= E_proc_B[b_min_id+k]-E_proc_B[b_min_id+k-1] ratio=val_delta/int_delta l_B11[j,1]=x[b_min_id+k]+ratio*(x[b_min_id+k]-x[b_min_id+k-1]) else: #No match, ignore l_B11[j,1]=np.nan # To find the crossing point consider l_A-l_B diff= l_A11[:,1]-l_B11[:,1] first_point=np.where( np.logical_not(np.isnan(diff)))[0][0] diff=diff[first_point:] cross=np.nan #Find the crossing point if first_point<0: if int(np.where(diff>0.0)[0].shape[0])>0: cross=int(np.where(diff<=0)[0][0]) else: if int(np.where(diff<0.0)[0].shape[0])>0: cross=np.where(diff>=0)[0][0] #Plot l_A-l_B ax=fig.add_subplot(311) ax.set_title(r"(a) Charging points in [{0:.0f},{1:.0f}] and discharging points in [{2:.0f},{3:.0f}]".format(x[a_min_id],x[a_max_id],x[b_min_id],x[b_max_id]),y=1.08) ax.plot(l_A11[:,0],l_A11[:,1]-l_B11[:,1], 'k-') ax.set_xlabel(r'Charge Demand, $u$ MW') ax.set_ylabel(r'$l_A(u)-l_B(u)$') if not np.isnan(cross): ax.plot(l_A11[first_point+cross,0],0.,'r+') ax.text(11000.,-80., r"$u=${0:.1f}, $l_A(u)=${1:.1f}, $l_B(u)=${2:.1f}".format(l_A11[first_point+cross,0],l_A11[first_point+cross,1],l_B11[first_point+cross,1])) ax.text(11000.,-85., r"$A$={0:,.1f}, $B$={1:,.1f}".format(F_proc_A[a_min_id+first_point+cross],F_proc_B[a_min_id+first_point+cross])) #Construct l_A and l_B for valid filling interval 2 and valid emptying interval 2 # Valid filling interval 2 [32571.0, 38357.0] index [12410,18196] # Valid emptying interval 2 [38358.0,54087.0] index [18197,33926] #This proves unsuccessful #Initialise arrays a_min_id=12410 a_max_id=18196 b_min_id=18197 b_max_id=33926 array_size = a_max_id-a_min_id l_A22=np.zeros((array_size,2)) l_B22=np.zeros((array_size,2)) #Run through possible charging points to construct l_A for j in range(array_size): i=a_min_id+j l_A22[j,0]=x[i] val_u=F_proc_A[i] #Does val_u lie in the range of possible values in the discharging region if val_u<=np.nanmax(E_proc_A[b_min_id:b_max_id]) and val_u>=np.nanmin(E_proc_A[b_min_id:b_max_id]): #find the value of l_A(u) k=np.where(E_proc_A[b_min_id:b_max_id]<=val_u)[0][0] val_delta=val_u-E_proc_A[b_min_id+k-1] int_delta= E_proc_A[b_min_id+k]-E_proc_A[b_min_id+k-1] ratio=val_delta/int_delta l_A22[j,1]=x[b_min_id+k-1]+ratio*(x[b_min_id+k]-x[b_min_id+k-1]) else: #No match, ignore l_A22[j,1]=np.nan #Run through possible charging points to construct l_B for j in range(array_size): i=a_min_id+j l_B22[j,0]=x[i] val_u=F_proc_B[i] #Does val_u lie in the range of possible values in the discharging region if val_u<=np.nanmax(E_proc_B[b_min_id:b_max_id]) and val_u>=np.nanmin(E_proc_B[b_min_id:b_max_id]): #find the value of l_B(u) k=np.where(E_proc_B[b_min_id:b_max_id]<=val_u)[0][0] val_delta=val_u-E_proc_B[b_min_id+k] int_delta= E_proc_B[b_min_id+k]-E_proc_B[b_min_id+k-1] ratio=val_delta/int_delta l_B22[j,1]=x[b_min_id+k]+ratio*(x[b_min_id+k]-x[b_min_id+k-1]) else: #No match, ignore l_B22[j,1]=np.nan # To find the crossing point consider l_A-l_B diff= l_A22[:,1]-l_B22[:,1] first_point=np.where( np.logical_not(np.isnan(diff)))[0][0] diff=diff[first_point:] cross=np.nan #Find the crossing point if first_point<0: if int(np.where(diff>0.0)[0].shape[0])>0: cross=int(np.where(diff<=0)[0][0]) else: if int(np.where(diff<0.0)[0].shape[0])>0: cross=np.where(diff>=0)[0][0] #Plot l_A-l_B ax=fig.add_subplot(312) ax.set_title(r"(b) Charging points in [{0:.0f},{1:.0f}] and discharging points in [{2:.0f},{3:.0f}]".format(x[a_min_id],x[a_max_id],x[b_min_id],x[b_max_id]),y=1.08) ax.plot(l_A22[:,0],l_A22[:,1]-l_B22[:,1], 'k-') ax.set_xlabel(r'Charge Demand, $u$ MW') ax.set_ylabel(r'$l_A(u)-l_B(u)$') if not np.isnan(cross): plt.plot(l_A22[first_point+cross,0],0.,'r+') plt.text(22000.,-80., r"$u=${0:.1f}, $l_A(u)=${1:.1f}, $l_B(u)=${2:.1f}".format(l_A22[first_point+cross,0],l_A22[first_point+cross,1],l_B22[first_point+cross,1])) plt.text(22000.,-85., r"$A$={0:,.1f}, $B$={1:,.1f}".format(F_proc_A[a_min_id+first_point+cross],F_proc_B[a_min_id+first_point+cross])) #Construct l_A and l_B for valid filling interval 1 and valid emptying interval 2 #Initialise arrays a_min_id=0 a_max_id=7962 b_min_id=18197 b_max_id=33926 array_size = a_max_id-a_min_id l_A12=np.zeros((array_size,2)) l_B12=np.zeros((array_size,2)) #Run through possible charging points to construct l_A for j in range(array_size): i=a_min_id+j l_A12[j,0]=x[i] val_u=F_proc_A[i] #Does val_u lie in the range of possible values in the discharging region if val_u<=np.nanmax(E_proc_A[b_min_id:b_max_id]) and val_u>=np.nanmin(E_proc_A[b_min_id:b_max_id]): #find the value of l_A(u) k=np.where(E_proc_A[b_min_id:b_max_id]<=val_u)[0][0] val_delta=val_u-E_proc_A[b_min_id+k-1] int_delta= E_proc_A[b_min_id+k]-E_proc_A[b_min_id+k-1] ratio=val_delta/int_delta l_A12[j,1]=x[b_min_id+k-1]+ratio*(x[b_min_id+k]-x[b_min_id+k-1]) else: #No match, ignore l_A12[j,1]=np.nan #Run through possible charging points to construct l_B for j in range(array_size): i=a_min_id+j l_B12[j,0]=x[i] val_u=F_proc_B[i] #Does val_u lie in the range of possible values in the discharging region if val_u<=np.nanmax(E_proc_B[b_min_id:b_max_id]) and val_u>=np.nanmin(E_proc_B[b_min_id:b_max_id]): #find the value of l_B(u) k=np.where(E_proc_B[b_min_id:b_max_id]<=val_u)[0][0] val_delta=val_u-E_proc_B[b_min_id+k] int_delta= E_proc_B[b_min_id+k]-E_proc_B[b_min_id+k-1] ratio=val_delta/int_delta l_B12[j,1]=x[b_min_id+k]+ratio*(x[b_min_id+k]-x[b_min_id+k-1]) else: #No match, ignore l_B12[j,1]=np.nan # To find the crossing point consider l_A-l_B diff= l_A12[:,1]-l_B12[:,1] first_point=np.where( np.logical_not(np.isnan(diff)))[0][0] diff=diff[first_point:] cross=np.nan #Find the crossing point if first_point<0: if int(np.where(diff>0.0)[0].shape[0])>0: cross=int(np.where(diff<=0)[0][0]) else: if int(np.where(diff<0.0)[0].shape[0])>0: cross=np.where(diff>=0)[0][0] #Plot l_A-l_B ax=fig.add_subplot(313) ax.set_title(r"(c) Charging points in [{0:.0f},{1:.0f}] and discharging points in [{2:.0f},{3:.0f}]".format(x[a_min_id],x[a_max_id],x[b_min_id],x[b_max_id]),y=1.08) ax.plot(l_A12[:,0],l_A12[:,1]-l_B12[:,1], 'k-') ax.set_xlabel(r'Charge Demand, $u$ MW') ax.set_ylabel(r'$l_A(u)-l_B(u)$') if not np.isnan(cross): #Crosses zero ratio= (0.0-diff[cross-1])/(diff[cross]-diff[cross-1]) a=l_A12[first_point+cross-1,0]+ratio*(l_A12[first_point+cross,0]-l_A12[first_point+cross-1,0]) b=l_A12[first_point+cross-1,1]+ratio*(l_A12[first_point+cross,1]-l_A12[first_point+cross-1,1]) A=F_proc_A[a_min_id+first_point+cross-1]+ratio*(F_proc_A[a_min_id+first_point+cross]-F_proc_A[a_min_id+first_point+cross-1]) B=F_proc_B[a_min_id+first_point+cross-1]+ratio*(F_proc_B[a_min_id+first_point+cross]-F_proc_B[a_min_id+first_point+cross-1]) ax.plot(a,0.,'k+', markersize=20,markeredgewidth=2) ax.text(a,-10., r"$a=${0:,.1f}, $b=${1:,.1f}".format(a,b)) ax.text(a,-20., r"$A$={0:,.1f}, $B$={1:,.1f}".format(A,B)) print("a ={0}, a_id={1}, b={2}, b_id={3}, A={4}, B={5}".format(a,a_min_id+first_point+cross-1,b,b_min_id+first_point+cross-1,A,B )) plt.tight_layout() fig.savefig('lAlB.pdf') plt.show() plt.close() np.save(str('l_A'),l_A12) np.save(str('l_B'),l_B12)