PyRAPS_SS.py

"""
09-2026

Copyright 2026, Battelle Energy Alliance, LLC, ALL RIGHTS RESERVED

@authors: Soumyadeep Nag, Shafiul Alam
"""

import csv
import psspy

#Desired head level
H0 = 0.9

#Input files
bsl_file = 'baseload.csv'
trate_file = 'gov_and_trate.csv'
derating_file = 'Affected_generators.csv'

Sid=-1   # All machines
Flag=1   # Only in-service machines at in-service plants
ierr, Nmach = psspy.amachcount(Sid, Flag)           # get no of machines in the subsystem
ierr, iMbus = psspy.amachint(Sid, Flag, 'NUMBER') # get machine bus numbers
ierr, cMids = psspy.amachchar(Sid, Flag, 'ID')    # get machine IDs   
iMbus=iMbus[0]
cMids=cMids[0]

#load flow parameters
ierr = psspy.solution_parameters_4(intgar2=200,realar5=.1)#,realar18 = .001

#read trate excel file--this files comes with the code
with open(trate_file, 'r') as file:
    csv_reader = csv.reader(file)
    header = next(csv_reader)  # Read the header row
    #print("Header:", header)
    rows_tr = list(csv_reader)
gt = [row[0] for row in rows_tr]
con_num = [row[1] for row in rows_tr]

#This function finds the maximum power for a given unit
def PM(ibus, genId,gt,con_num):
    ierr,Pm = psspy.macdat(ibus, genId, 'PMAX')
    ierr,Mb = psspy.macdat(ibus, genId, 'MBASE')
    ierr, gov_mdl = psspy.mdlnam(ibus, genId, 'GOV')
    if not gov_mdl is None:
        gov_mdl = gov_mdl.strip()
    if gov_mdl in gt:
        ierr, icon0 = psspy.mdlind(ibus, genId, 'GOV', 'CON') # get initial CON address (index)
        ierr,Trate=psspy.dsrval('CON', icon0+int(con_num[gt.index(gov_mdl)])-1)
    else:
        Trate = 1000000000000000
    
    PMAX = min([Trate,Pm,Mb])
    
    return PMAX
    
with open(bsl_file, 'r') as file:
    csv_reader = csv.reader(file)
    header = next(csv_reader)  # Read the header row
    #print("Header:", header)
    rows_bsl = list(csv_reader)

bus_num = [row[0] for row in rows_bsl]
bus_id = [row[2] for row in rows_bsl]
bsl_flag = [row[3] for row in rows_bsl]
  
with open(derating_file, 'r') as file:
    csv_reader = csv.reader(file)
    header = next(csv_reader)  # Read the header row
    #print("Header:", header)
    rows = list(csv_reader)
dert_bus = [row[1] for row in rows]

#This function performs redispatch 
def adjustgen(Nmach,iMbus,cMids,Pdiff,bus_num,bus_id,bsl_flag):
    
    with open ('reserve_report.csv', mode='w', newline='') as file:
        writer = csv.writer(file)
        # Write header
        header = ['bus name']+['Gen ID on the bus']+['Reserve_up']+['Reserve_dn']+['Gen_old']+['Gen_new']+['Pmax']+['Pmin']
        writer.writerow(header)
        sys_res_up = 0
        sys_res_dn = 0
        #find system reserve
        for i in range(0,len(bus_num)):
            for iM in range(0,Nmach): 
                ibus=iMbus[iM]
                genId=cMids[iM]
                genId = genId.strip()
                ierr, status = psspy.macint(ibus,genId,'STATUS')
                ierr, genMdl = psspy.mdlnam(ibus,genId,'GEN')
                if str(ibus) == bus_num[i] and genId == bus_id[i] and (status == 1) and (not genMdl is None) and bsl_flag[i]=='0':#in bus_num:    
                    ierr,Pg = psspy.macdat(ibus, genId, 'P')
                    Pmax = PM(ibus, genId,gt,con_num)
                    ierr,Pmin = psspy.macdat(ibus, genId, 'PMIN')
                    if Pg>0:
                        if (Pmax-Pg)>10:
                            sys_res_up += Pmax-Pg
                            sys_res_dn += Pg-Pmin
                        
                    break
        
        sys_res_up_new = 0
#        if -Pdiff>sys_res_up:
#            Pdiff = -.9*sys_res_up
#            print('reserve not sufficient')
        
        #adjust generation
        for i in range(0,len(bus_num)):
            for iM in range(0,Nmach): 
                ibus=iMbus[iM]
                genId=cMids[iM]
                genId = genId.strip()
                ierr, status = psspy.macint(ibus,genId,'STATUS')
                ierr, genMdl = psspy.mdlnam(ibus,genId,'GEN')
                if str(ibus) == bus_num[i] and genId == bus_id[i] and (status == 1) and (not genMdl is None) and int(bsl_flag[i])==0:#in bus_num:  
                    ierr,Pg_old = psspy.macdat(ibus, genId, 'P')
                    Pmax = PM(ibus, genId,gt,con_num)
                    ierr,Pmin = psspy.macdat(ibus, genId, 'PMIN')
                    if Pg_old>0:
                        if Pdiff<0:
                            if (Pmax-Pg_old)>10:
                                reserve = Pmax-Pg_old
                                Pg_new = Pg_old - ((reserve/sys_res_up)*Pdiff) #reserve of machine in consideration
                                if Pg_new>=Pmax:
                                    Pg_new = Pmax
                                row = [ibus]+[genId]+[Pmax-Pg_old]+[Pg_old-Pmin]+[Pg_old]+[Pg_new]+[Pmax]+[Pmin]
                                ierr = psspy.machine_data_4(ibus,genId,realar1=Pg_new)
                                sys_res_up_new += Pmax-Pg_new
                                writer.writerow(row)
                                break
                            else:
                                Pg_new = Pg_old
                        elif Pdiff>0:
                            reserve = Pg_old-Pmin
                            Pg_new = Pg_old - ((reserve/sys_res_dn)*Pdiff) #reserve of machine in consideration
                            row = [ibus]+[genId]+[Pmax-Pg_old]+[Pg_old-Pmin]+[Pg_old]+[Pg_new]+[Pmax]+[Pmin]
                            ierr = psspy.machine_data_4(ibus,genId,realar1=Pg_new)
                            writer.writerow(row)
                            break

    return sys_res_up_new

#adjust Pmax --- This is not used 
def adjustpmax(Nmach,iMbus,cMids,Pmax_diff,bus_num,bus_id,bsl_flag,dert_bus,Pgen_diff):
    
    sys_res_up = 0
    sys_res_dn = 0
    #find system reserve
    for iM in range(0,Nmach): 
        ibus=iMbus[iM]
        genId=cMids[iM]
        if str(ibus) in bus_num:
            ind = bus_num.index(str(ibus))
            genId = genId.strip()
            if genId == bus_id[ind]:
                if int(bsl_flag[ind])==0:
                    ierr,Pg = psspy.macdat(ibus, genId, 'P')
                    Pmax = PM(ibus, genId,gt,con_num)
                    ierr,Pmin = psspy.macdat(ibus, genId, 'PMIN')
                    sys_res_up += Pmax-Pg
                    sys_res_dn += Pg-Pmin
                        
    if sys_res_up<abs(Pgen_diff):
        Pmax_diff = -(abs(Pgen_diff)-sys_res_up)
        sys_max = 0
        #find system reserve
        for iM in range(0,Nmach): 
            ibus=iMbus[iM]
            genId=cMids[iM]
            if str(ibus) in bus_num:
                if not str(ibus) in dert_bus:
                    ind = bus_num.index(str(ibus))
                    genId = genId.strip()
                    if genId == bus_id[ind]:
                        if int(bsl_flag[ind])==0:
                            Pmax = PM(ibus, genId,gt,con_num)
                            sys_max += Pmax
                        
        #adjust Pmax
        for iM in range(0,Nmach): 
            ibus=iMbus[iM]
            genId=cMids[iM]
            if str(ibus) in bus_num:
                ind = bus_num.index(str(ibus))
                if not str(ibus) in dert_bus:
                    genId = genId.strip()
                    if genId == bus_id[ind]:
                        if int(bsl_flag[ind])==0:
                            ierr,Pmax_old = psspy.macdat(ibus, genId, 'PMAX')
                            Pmax_new = Pmax_old - ((Pmax_old/sys_max)*Pmax_diff) #reserve of machine in consideration
                            ierr = psspy.machine_data_4(ibus,genId,realar5=Pmax_new)
                            print('Pmax changed')
                        
Pgen_diff = 0
Pmax_diff = 0

#Perform derating and Calculate power deficit
for i in range(len(rows)):
   # i = 20
    x = rows[i]
    bus_number = int(x[1])
    gen_id = x[3]

    print("BUS_NUMBER:", bus_number)
    print("GEN_ID:", gen_id)
   
    #change Pgen and Pmax new of hydro unit
    ierr, status = psspy.macint(bus_number, gen_id, 'STATUS')   
    if status == 1:
        
        ierr,Pgen_old = psspy.macdat(bus_number, gen_id, 'P')
        ierr,Pmin_old = psspy.macdat(bus_number, gen_id, 'PMIN')
        Pmax_old = PM(bus_number, gen_id,gt,con_num)
        
        if Pgen_old>0:
            
            ierr, gov_mdl = psspy.mdlnam(bus_number, gen_id, 'GOV') # get governor model name
            
            if not gov_mdl is None:
                gov_mdl = gov_mdl.strip()
            
            if gov_mdl == 'HYG3U1':
                ierr, icon0 = psspy.mdlind(bus_number, gen_id, 'GOV', 'CON') # get initial CON address (index)
                ierr,g_max=psspy.dsrval('CON', icon0+13-1)
                ierr,qnl=psspy.dsrval('CON', icon0+29-1)
                ierr,At=psspy.dsrval('CON', icon0+31-1)
            
                #Pmax_new = Pmax_old*H0**(3/2)
                Pmax_new = Pmax_old*At*((g_max*(H0**1.5))-(qnl*H0))
                Pmin_new = Pmin_old
                Pgen_new = Pgen_old
            else:
                Pmax_new = Pmax_old*H0**(3/2)
                Pmin_new = Pmin_old
                Pgen_new = Pgen_old
            
            if Pmax_new<Pmin_new:
                Pmax_new = Pmin_new
                Pgen_new = Pmax_new
                
            if Pgen_new>Pmax_new:
                #Pmax_new = Pgen_new #POM version
                Pgen_new = .95*Pmax_new#leave space for losses and response

            if Pgen_new<Pmin_new:
                Pgen_new = Pmin_new

            #change hydro plant generation 
            ierr = psspy.machine_data_4(bus_number, gen_id,realar1=Pgen_new,realar5=Pmax_new,realar6=Pmin_new)
            
            #find difference in generation
            Pgen_diff += Pgen_new-Pgen_old
            Pmax_diff += Pmax_new-Pmax_old
    else:
        print('Machine in deration list is out of service')  
   

#adjust other non-renewable and available generators 
SUR = adjustgen(Nmach,iMbus,cMids,Pgen_diff,bus_num,bus_id,bsl_flag)

                    
# Solve the load flow
psspy.fnsl([0,0,0,0,0,0,0,0])  # Full-Newton load flow solution

# Check for convergence
solved1 = psspy.solved()

#redispatch slack
if solved1 == 0:
    print("First Load flow solved successfully. Adjusting swing bus generation.")
    print('Total system upward reserve:',SUR)
    print('Total adjustment needed for Pmax:',Pmax_diff)
    print('Total adjustment needed for Pgen:',Pgen_diff)
    
    for iM in range(0,Nmach): 
        ibus=iMbus[iM]
        genId=cMids[iM]
        ierr, ival = psspy.busint(ibus,'TYPE')
        if ival==3:
            ierr,Pgen_old = psspy.macdat(ibus, genId, 'P')
            ierr,Pmin_old = psspy.macdat(ibus, genId, 'PMIN')
            Pmax_old = PM(ibus, genId, gt, con_num)
            SW_bus = ibus
            if Pgen_old>Pmax_old:
                SW_bus = ibus
                Pgen_diff = 1.2*(Pmax_old-Pgen_old)
                Pmax_diff = 0
                    
                #adjust other non-renewable and available generators 
                #adjustpmax(Nmach,iMbus,cMids,Pmax_diff,bus_num,bus_id,bsl_flag,dert_bus,Pgen_diff)

                #adjust other non-renewable and available generators 
                SUR = adjustgen(Nmach,iMbus,cMids,Pgen_diff,bus_num,bus_id,bsl_flag)
            else:
                print("Swing bus adjustment was not required")


    # Solve the load flow
    psspy.fnsl([0,0,0,0,0,0,0,0])  # Full-Newton load flow solution

    # Check for convergence
    solved2 = psspy.solved()

#save new case files
if solved1 == 0:
    print("First load flow converged successfully.")
    if solved2==0:
        print("Second Load flow solved successfully.")
        
        file = "HS_NW_CAL90.sav"
        psspy.save(file)
    else:
        print("Second Load flow did not converge.")
else:
    print("First Load flow did not converge.")

print('Total system upward reserve:',SUR)
print('Total adjustment needed for Pmax:',Pmax_diff)
print('Total adjustment needed for Pgen:',Pgen_diff)
print('Swing bus is:',SW_bus)
print("Modified case is save in: ", file)