fix bug

src/necost/__init__.py

@@ -51,7 +51,7 @@ def __init__(
         self.pv_internals_ref = pv_internals_ref
         self.pv_turb_ref = pv_turb_ref
         self.m_turb = m_turb
-        self.number_rods_ref = assembly_ref * 264,
+        self.number_rods_ref = assembly_ref * 264
 
         # Modifications to the values in the Monte Carlo samples
         self.core_power = self.data["core_pwr_func_ref"] * self.data["ref_therm_pwr"]
@@ -83,21 +83,19 @@ def run(self):
                 self.data["fuel_density"] * self.data["HM_weight_percent"]),
             inplace=True
         )
-
         # reference core specific power W/gHM.
-        q_sp_ref = self.data["ref_therm_pwr"] / self.number_rods_ref * self.l_h_ref * np.pi * (
+        mass_hm_core_ref = self.number_rods_ref * self.l_h_ref * np.pi * (
             self.fuel_diam_ref / 2) ** 2 * (self.fuel_density_ref * self.weight_hm_ref)
-
+        q_sp_ref = self.data["ref_therm_pwr"] / mass_hm_core_ref
         q_sp = self.core_power / self.data["HM_mass_direct_spec"]  # based on whole core -  W/gHM.
         efpy_ref = (self.ref_burnup * (1000 / 365.25)) / q_sp_ref  # for reference core [y]
         efpy = (self.data["max_burnup"] * (1000 / 365.25)) / q_sp  # EFPY at total discharge burnup [years]
 
         # --------------------------------  Interest rate-related parameters
         lev_factor = self.data["discount_rate"] / (1 - np.exp(-self.data["discount_rate"] * self.t_plant))
-        e_year = self.core_power * self.data["therm_efficiency"] * (self.data["L_direct_spec"] / 1000) * 8766
+        e_year = self.core_power * self.data["therm_efficiency"]/100 * (self.data["L_direct_spec"] / 1000) * 8766
 
         l_inp, n_cycl, t_cyc, planned_outage = function_1(self.data, self.batches_ref, self.t_plant, efpy, efpy_ref)
-
         m_fabrication, m_enrichment, pv_front_end_u_1 = front_end_1(self.data)
         pv_front_end_r_1 = front_end_2(self.data)
         levelized_front_end = scale_up(
@@ -110,10 +108,8 @@ def run(self):
             n_cycl,
             t_cyc
         )
-
         levelized_back_end = backend(self.data, self.t_plant, lev_factor, e_year, n_cycl, t_cyc)
         levelized_fcc = levelized_front_end + levelized_back_end
-
         levelized_hydride_om = levelized_om_cost(
             self.data,
             self.t_plant,
@@ -124,7 +120,6 @@ def run(self):
             t_cyc,
             planned_outage
         )
-
         # --------------------------------  capital costs
         levelized_hydride_cap, levelized_hydride_cost = capital_cost(
             self.data,

src/necost/backend.py

@@ -51,26 +51,22 @@ def backend(data: pd.DataFrame, t_plant, lev_factor, e_year, n_cycl: np.ndarray,
 
     # Most create a copy of the values. Otherwise, `pv_back_end_o` will also be updated below.
     pv_back_end_o = pv_back_end.copy(deep=True)
-
     # TODO: Should the `np.arange(...)` start at 1?
-    k = n_cycl - 1
-    indices = np.arange(np.max(k))  # Generate an array of indices from 0 to MAX(n_cycl - 2)
-    m = np.where(indices < k[:, np.newaxis], indices, -np.inf)  # Apply the mask to fill the matrix.
-
+    # k = n_cycl - 1
+    indices = np.arange(1, np.max(n_cycl))  # Generate an array of indices from 0 to MAX(n_cycl - 2)
+    m = np.where(indices < n_cycl[:, np.newaxis], indices, -np.inf)  # Apply the mask to fill the matrix.
     # We used -np.inf to define `matrix` above so that those values become zero after exponentiating,
     # and so the sum will not be affected by those values. If we used 0 instead of -inf, those values
     # would become 1 after exponentiation, which would be added in the sum.
     pv_back_end *= 1 + np.exp(
         (data["escalation_rate_back"] - data["discount_rate"]).values.reshape(-1, 1) * m * t_cyc.reshape(-1, 1)
     ).sum(axis=1)
-
     # as above, but this is the residual of the last cycle, in percent over the cycle (~ 15 #) # Not realistic
     # because of transportation costs (James of Exelon): only the fraction of t_cyc utilized is a cost,
     # while the rest is sold.
     pv_back_end += (pv_back_end_o * np.exp(
         (data["escalation_rate_back"] - data["discount_rate"]) * n_cycl * t_cyc
     )) * (t_plant - (n_cycl * t_cyc)) / t_cyc
-
     # adds the 1 mills/kWh or not depending if it has been requested
     # (on what fraction??, in this case total electricity produced, not core discharged)
     levelized_back_end = np.where(

src/necost/costs.py

@@ -30,24 +30,20 @@ def levelized_om_cost(
     pv_om = np.where(
         data["OM_fixed_cost"] <= 0,
         cf_ro + cf_fo + cf_pers,
-        data["OM_direct_spec"] * (core_power / 1000) * data["therm_efficiency"]
+        data["OM_direct_spec"] * (core_power / 1000) * data["therm_efficiency"]/100
     )
-
     # PV of O&M costs at t = 1 year ($).
     # TODO: Should the `np.arange(...)` start at 1?
-    q, _ = np.meshgrid(np.arange(t_plant - 1), range(len(data)))
+    q, _ = np.meshgrid(np.arange(1, t_plant ), range(len(data)))
     epsilon = np.exp((data["escalation_rate_OM"] - data["discount_rate"]).values.reshape(-1, 1) * q)
     pv_om = pv_om * (1 + epsilon.sum(axis=1))
-
     # Total PV of O & M costs at t = 0 year ($)
     pv_om = pv_om * np.exp(-data["discount_rate"] * 1)
-
     # levelized O & M cost over life of plant ($/year)
     levelized_om_total = pv_om * lev_factor
 
     # levelized O & M unit cost over life of plant (mills/kW-hre)
     levelized_hydride_om = levelized_om_total * 1000 / e_year + o_m_fixed
-
     return levelized_hydride_om
 
 
@@ -71,15 +67,15 @@ def compute_pv_cap(
         constrct_years = data["constrct_years"].values
         capital_cost_vals = data["capital_cost"].values
         ref_therm_pwr = data["ref_therm_pwr"].values
-        therm_efficiency = data["therm_efficiency"].values
+        therm_efficiency = data["therm_efficiency"].values/100
         escalation_cost_pwr = data["escalation_cost_pwr"].values
         apr_r_constr = (1 + interest_rate_constrct / 4) ** 4 - 1
 
         # To compute the `q_cost_no_interest` outside the loop, we must use this masking
         # technique since `constrct_years` may have different values, leading to different
         # lengths for the `q_cost_no_interest` matrix.
-        indices = np.arange(0, np.max(constrct_years), 0.25)  # Generate an array of indices
-        m = np.where(indices < constrct_years[:, np.newaxis], indices, np.nan)  # Apply the mask to fill the matrix
+        indices = np.arange(0, np.max(constrct_years+0.25), 0.25)  # Generate an array of indices
+        m = np.where(indices < constrct_years[:, np.newaxis]+0.25, indices, np.nan)  # Apply the mask to fill the matrix
         x_constr = -np.pi / 2 + np.pi * (m / constrct_years.reshape(-1, 1))
         s_curve_cum = 0.5 * (np.sin(x_constr) + 1)
         s_curve_diff = np.diff(s_curve_cum)
@@ -90,11 +86,10 @@ def compute_pv_cap(
         tot_capital_cost = np.empty(len(data))  # --> populated by the loop
         for idx in range(len(interest_rate_constrct)):
             # Select the correct values for the current computation. Values outside the slice are NaN.
-            current_q_cost_no_interest = q_cost_no_interest[idx, :(int(constrct_years[idx] / 0.25) - 1)]
+            current_q_cost_no_interest = q_cost_no_interest[idx, :(int(constrct_years[idx] / 0.25) )]
 
             tot_capital_cost[idx] = current_q_cost_no_interest @ np.full(len(current_q_cost_no_interest), (
-                1 + apr_r_constr[idx] / 4)) ** (np.arange(constrct_years[idx] * 4, 1, -1) - 0.5)
-
+                1 + apr_r_constr[idx] / 4)) ** (np.arange(constrct_years[idx] * 4, 0, -1) - 0.5)
         # If cap cost exp is 1 ==> ref_therm_pwr cancels out and only core_power is used
         pv_cap = tot_capital_cost * ref_therm_pwr / 1000 * therm_efficiency * (
             core_power / ref_therm_pwr) ** escalation_cost_pwr

src/necost/frontend.py

@@ -166,22 +166,19 @@ def scale_up(
     pv_front_end_tot = pv_front_end_u_1 * data["frac_core_loaded_nat"] + pv_front_end_r_1 * data[
         "frac_core_loaded_reprcsd"]
     pv_front_end = pv_front_end_tot * (data["HM_mass_direct_spec"] / 1000 / data["num_batches"])
-
     #  present value of the entire front end fuel cycle costs for all batches throughout the core lifetime
     pv_front_end_o = pv_front_end.copy(deep=True)
 
     # TODO: Should the `np.arange(...)` start at 1?
-    k = n_cycl - 1
-    indices = np.arange(np.max(k))  # Generate an array of indices from 0 to N-1
-    m = np.where(indices < k[:, np.newaxis], indices, 0)  # Apply the mask to fill the matrix
-
+    # k = n_cycl - 1
+    indices = np.arange(1,np.max(n_cycl))  # Generate an array of indices from 0 to N-1
+    m = np.where(indices < n_cycl[:, np.newaxis], indices, 0)  # Apply the mask to fill the matrix
     # We used -np.inf to define `matrix` above so that those values become zero after exponentiating,
     # and so the sum will not be affected by those values. If we used 0 instead of -inf, those values
     # would become 1 after exponentiation, which would be added in the sum.
     pv_front_end *= (1 + np.exp(
         (data["escalation_rate_front"] - data["discount_rate"]).values.reshape(-1, 1) * m * t_cyc.reshape(-1, 1)
     ).sum(axis=1))
-
     # as above, but this is the residual of the last cycle, in percent over the cycle (~ 15 #)
     # Not realistic because of transportation costs (James of Exelon): only the fraction of t_cyc
     # utilized is a cost, while the rest is sold.

tutorial/cal_lcae.py

@@ -56,7 +56,10 @@
 # result_r3 = necost_r3.run()
 result_r5 = necost_r5.run()
 # result_r10 = necost_r10.run()
-print(result_r5)
+# transpose result data
+result_t = result_r5.transpose()
+print(result_t)
+