A Differentiable Model for Optimizing Hybridization of Industrial Process Heat Systems with Concentrating Solar Thermal Power

A dynamic model of a concentrating solar thermal array and thermal energy storage system is presented that is differentiable in the design decision variables: solar aperture area and thermal energy storage capacity. The model takes as input the geographic location of the system of interest and the corresponding discrete hourly solar insolation data, and calculates the annual thermal and economic performance of a particular design. The model is formulated for use in determining optimal hybridization strategies for industrial process heat applications using deterministic gradient-based optimization algorithms. Both convex and nonconvex problem formulations are presented. To demonstrate the practicability of the models, they were applied to four different case studies for three disparate geographic locations in the US. The corresponding optimal design problems were solved to global optimality using deterministic gradient-based optimization algorithms. The model and optimization-based analysis provide a rigorous quantitative design and investment decision-making framework for engineering design and project investment workflows.


Introduction
Industrial process systems are enormous consumers of energy (9.23 TWh th ), accounting for roughly one-third of all delivered energy use (28.64 TWh th ) in the United States in 2017 [1].Roughly 60% of the energy consumption in the industrial sector is from the manufacturing sector, which accounts for roughly 20% of the total US energy consumption [2].The industrial sector has historically been the largest consumer of natural gas [1] as both a feedstock and a primary energy source.In manufacturing, industrial process heat (IPH) accounts for roughly 70% of the total process energy consumed in the US [3].When broken down by fuel consumption, the manufacturing sector accounts for roughly 95% of all IPH demand in the US [3].Further, since the majority of these process systems require low-or medium-temperature process heat (<250 • C) [3,4], there is a substantial opportunity to reduce fossil fuel consumption and carbon emissions by augmenting existing IPH systems with solar thermal power and designing new systems incorporating solar thermal power.
Numerous studies have been conducted over the past 40 years on the potential that solar energy can play within IPH systems [4][5][6][7][8][9].Active interest in the topic in recent years has been motivated largely by the concern over unsustainable energy consumption, greenhouse gas emissions, and overall environmental impact of commercial and industrial sectors [10].Attention to the feasibility of solar applications to specific industries has been made, such as: food and beverage [11], mining [12,13], agriculture [14,15], textiles [11,[16][17][18][19], pulp and paper [20,21], chemicals [22,23], and building construction materials [4], among others.Recent reviews [4,6,7,9] resolve the individual industries by their similar unit operations (e.g., drying, hot water, sterilization, etc.).Results of such studies indicate general feasibility of solar in IPH applications; however, the economic viability is extremely sensitive to the application and geographic region of the project.Given the diversity of IPH demand profiles across industries and applications as well as the widely-varying solar resource performance profiles across geographic locations, a formal numerical simulation and optimization-based design approach is necessary to ensure the economic viability of hybrid solar-IPH systems.
A formal mathematical optimization approach was applied by Ghobeity and Mitsos [28] to the design and operation of a solar thermal receiver and thermal energy storage (TES) system with a constant IPH demand for high-temperature steam generation (for power cycles).A sequential design approach was applied that first considers minimizing total solar cost with the solar array size and TES capacity as design decision variables and several operating decision variables related to the proposed dynamical model [28].The model was developed for use with gradient-based local optimization solvers and utilized a total daily concentrated energy value based on a single winter day and a heliostat field efficiency calculated for a location in Cyprus [28].The economic objective was considered as a simple linear cost model with fixed specific costs for the TES medium and solar array.The second optimization problem considered the yearly operation for the fixed design based on the winter day.Thus, for the fixed design, feasible optimal operating conditions were found with respect to 12 representative days throughout the year [28].
Modeling and simulation for hybridizing power cycles with concentrating solar thermal (CST) power was recently considered by Gunasekaran and coworkers [30] with system optimization by Ghasemi and coworkers [29] and Ayub and coworkers [36].In [29], the focus was on optimizing the operation of the hybrid system without considering TES.In [36], economics were considered to optimize the levelized cost of electricity.However, due to the regional dependence of the solar model and economic assumptions, optimization with respect to the economics of the parabolic trough solar concentrator (PTC) array sizing were omitted.Additionally, in [36], no TES was considered.The thermo-economic model was implemented in Engineering Equation Solver (EES) and a shortcut method is used to optimize the hybrid system using intermediate variables calculated from the AspenPlus sequential-quadratic programming (SQP) optimizer.
A simplified approach to the optimal design of solar water heating systems was presented by Yan et al. [33].The authors developed a high-level energy conservation model of the water heating system including TES as a hot water tank.The design decision variables were taken as the array size and the TES capacity and the objective was the lifecycle savings.The authors [33] proposed an approach based on diminishing marginal utility and identified optimal solutions based on the intersection of marginal energy savings and marginal embodied energy curves.Although their simplified approach was effective for the water-heating application, IPH applications and their dramatically differing energy demand profiles complicate the optimal design problem as the total utilizable energy and lifecycle savings exhibit large degrees of nonlinear sensitivity to the decision variables [33].
A hybrid solar-IPH design strategy based on heat integration was considered by Baniassadi et al. [34].The authors considered fixed array designs determined by maximizing the solar fraction for a particular pre-determined investment and analyzed the various options for coupling solar energy with process streams.They proposed an algorithm which uses the various process stream temperatures and pinch analysis to identify the best process stream to exchange solar energy with.They utilized a commercial simulator to calculate the performance of the solar system and subsequently the economics of the retrofit strategy [34].Although the analysis identifies some technical challenges of physically coupling the solar thermal system with the process from a heat integration perspective, it does not adequately address the fundamental design considerations of the solar system itself [34].Further, their approach cannot be used within mathematical optimization-based design strategies as the solar performance model is only accessible as an ad-hoc black-box simulation [34].
In [31], Silva and coworkers considered a formal simulation and optimization-based design approach to hybrid solar-IPH systems.In their work, a single-tank thermocline TES device was modeled along with a PTC array for medium-temperature (100-250 • C) IPH applications.Three economic objective functions were considered: levelized cost of energy, payback time, and lifecycle savings as their base case [31].Four design variables were considered in their work: the number of PTC modules in series, the number of PTC modules in parallel, the spacing between the rows, and the TES capacity [31].The authors propose a derivative-free non-deterministic memetic optimization algorithm based on a modified genetic algorithm [31].They make explicitly clear that this approach is in response to the major challenges in optimizing dynamical systems based on simulation: possible discontinuities, nonlinearity, and expensive complex objective function evaluations [31].
In [32], a hybrid solar-IPH system for zero-liquid discharge desalination (<250 • C) was considered.Similar to [31], a formal simulation and optimization-based design approach was taken using the net-present value as the economic objective.The design variables considered for the base-case were the number of PTC modules and the TES capacity.A formal robust design problem was considered as a nonconvex non-cooperative game by taking into account uncertainty in the natural gas and product markets and solving the problem formulated as a semi-infinite program [32].Again, the authors note that due to the non-differentiable nature of the model, a derivative-free genetic algorithm-based optimization approach must be used with subsequent analysis to verify feasibility and optimality of solutions [32].
The past approaches to modeling, simulation, and optimization of solar-IPH systems have all shown favorable results using model formulations or approaches specific to the application of interest or more general formulations which pose challenges for deterministic gradient-based optimization.There exists a need for a dynamic model of a CST-TES system for designing optimal hybrid solar-IPH systems with accurate economic models for assessing technology investment risk and the economic viability of such projects.In this paper, a model is developed for low-to medium-temperature CST-TES systems for assessing the economic viability and overall project feasibility of hybrid solar-IPH systems.The model accounts for a single-axis tracking PTC array, a TES device, and utilizes high-accuracy hourly solar resource data for the United States available from an open-access database.The model is developed specifically for use with deterministic (global) optimization algorithms, and therefore accounts for convex analysis, differentiability, and the trade-offs between computational efficiency and model accuracy.
In Section 2, the hybrid solar-IPH process performance model will be developed which includes a smooth formulation, a model of a PTC with solar angle calculations for high-resolution dynamic simulation of CST performance, an economic model for assessing design feasibility and investment decision-making, and a mathematical optimization formulation.Additionally, the key theoretical results and analysis of the model are presented in Section 2. In Section 3, four numerical case studies for varying fuel costs and IPH demand profiles across three geographically-disparate locations are presented and discussed followed by the Conclusion in Section 4.

Concentrating Solar Thermal Model
For a surface or defined geometry oriented horizontally, the performance of the solar energy collector relies on the surface tilt angle and the incidence angle of the irradiation striking the surface.Here, we consider the model of a PTC, illustrated in Figure 1 with the aperture area of a single trough given by A a = L a L m .The total optical efficiency of the reflector is defined in the following.Definition 1 (Optical Efficiency, η opt [30]).The optical efficiency of the PTC η opt is given by: where i ∈ [0, 1] for i = 1, . . ., 6 are efficiency factors related to various components of the PTC and the incidence angle modifier K is given by the empirical relationship [39,40]: Here, 1 is the shadowing factor of the receiver tube (=0 for fully-shadowed and =1 for no shadowing), 2 is the tracking error factor (=1 for no tracking error), 3 is the geometry error (=1 for no geometric or mirror alignment errors), 4 is the mirror dirt factor (=1 for perfectly clean mirrors), 5 is the dirt factor for the glass envelope enclosing the receiver tube (=1 for perfectly clean glass), and 6 are the unaccounted for losses (=1 for no unaccounted losses).The factor r m is the reflectance of perfectly clean mirrors (=1 for 100% reflectance), a r is the absorptance of the receiver tube itself (=1 for a perfect black-body), and τ is the transmittance of the glass envelope enclosing the receiver tube.The factor K accounts for the fact that irradiance is not normal to the PTC aperture (unless the array has two-axis tracking).The incidence angle θ follows from the solar position model in Appendix A.
With the above quantities, the specific thermal power potential (in units of kW th /m 2 ) is given by q = I d η opt /1000, with I d the direct-normal irradiance (DNI) at a particular moment in time and geographic location (in standard units of W/m 2 ).Note the DNI values I d come from experimental measurements.These values are available for specific geographic locations in the US (and some international regions) for annual time horizons with hourly resolution (i.e., 8760 discrete hourly values), via the National Renewable Energy Laboratory's (NREL) National Solar Radiation Database (NSRDB) [41].In this work, a geographic location is specified and the typical meteorological year (TMY) data is downloaded.The numerical values of the PTC modeling parameters used in this study can be found in Table A1 of Appendix B. Finally, the solar thermal power produced (in units of kW) by a solar array is given by qs = qx a , with x a ∈ R + as the total aperture area of the solar array.

Process System Model
The process system of interest is represented by the block-flow diagram in Figure 2. Physically, the system consists of the solar energy collector as the source of solar thermal energy for the system, the TES, the IPH system which exists as an energy demand, and the conventional (fossil) energy source for times of inadequate solar energy.The process system model is developed from a high-level with the following simplifying assumptions: 1. Heat losses to the environment from piping and heat exchangers is negligible.2. Heat losses to the environment from the TES is negligible.3. The process system requires low to medium temperature IPH (i.e., ≤250 • C). 4. Heat is always available at or above the minimum temperature as required by the IPH system. 5.The IPH system energy demand is not dependent on the state or design decisions of the solar energy system.
Under these assumptions, the process system model is developed from the overall energy balance on the entire system: qs where qs is the instantaneous thermal power provided by the solar system, qng is the instantaneous thermal power provided by the conventional fossil-fuel (e.g., natural gas) heating system, qp is the instantaneous thermal power demand by the industrial process (IPH), qts is the instantaneous power supply/demand of the solar thermal energy storage system, and ql is the lost solar thermal power due to meeting the IPH requirements and reaching the maximum thermal storage capacity (i.e., the unutilizable solar energy).Note that qts can be both positive (charging state) and negative (discharging state).The complicating details for accounting for energy in the system comes from the capacity limitations of the TES.For example, without TES, the solar array will simply deliver all its energy to the IPH system until the full power demand qp is met.At which point, the solar array will simply de-focus, or we account for the positive difference between the supply and demand as losses or unutilizable solar ql .The negative difference between the supply and demand is simply accounted for by the conventional heating system qng .However, with the existence of a TES device, the positive difference between solar thermal power supply and IPH demand can be stored and used in future times when the solar system alone cannot meet the IPH demand.This complicates the model since we assert a preference for solar energy over conventional energy sources, and hence must account for the finite capacity of the storage device and its stored capacity at each instant in time.
Rearranging Equation ( 4) to focus on the storage device, we get: qts + qa = qs − qp where we have introduced qa = ql − qng for ancillary energy accounting.Physically, this quantity is used for load balancing control (i.e., balancing energy supply and demand).Equivalently, we write: assuming the TES device is empty at the start of the simulation.Applying the explicit Euler integration scheme with h = ∆t as the discrete time step-size (1 h for standard solar data), we can write: Next, we develop the model to account for the finite capacity of the TES device.Since we know that the energy capacity of the TES device at any instant can only take values between zero and the maximum capacity, we write: where x ts is the storage capacity variable scaled in hours of peak process demand qpeak p = max i qi p , determined by the process operations schedule model (known a priori).The inner max binary operation accounts for the condition that if there was excess energy available from the combined solar array and TES at the previous time, then we should consider storing it for future use.The outer min binary operation accounts for the condition that, if we are considering storing energy, we can only store at or below the maximum TES capacity.Equivalently, we can write min qpeak p x ts , max{0, where the mid operator selects the median value of the three arguments.Combining Equations ( 5)-( 7), the overall energy balance becomes: From this model, there are three possible states from the load-balancing perspective at each discrete time point i: The PTC is not providing enough energy to meet the demand of the process (TES moves into the discharging state) and there is not enough energy in the TES to meet the demand of the process over this time step.In this case, conventional energy must be supplied to the process to meet the full demand of the process.

( qi
The PTC is providing enough energy to meet the demand of the process (TES moves into the charging state) and there is not enough available capacity to store all excess energy in the TES device.In this case, the excess energy must be rejected.

( qi
The TES is either in the discharging state with enough stored energy to meet the full demand, or it is in the charging state and has enough available capacity to store the solar energy in excess of the process demand.
With these relationships, we can write: so that ql and qng are accurately accounted for.At any given discrete time point, we have the following relationship for load balancing: Note that losses are accounted for here as simply the unutilizable solar power at any given time step and do not account for standard thermal losses to the environment which occur in thermal systems.As stated previously, heat losses for the TES are expected to be minimal over the timescale of charging/discharging cycles.This is valid for well-insulated systems that undergo complete or nearly-complete discharge over the course of the day.For systems that have long-term (i.e., many days or months) storage, this assumption will not be valid.Further, TES systems to be installed outdoors and above-ground may need special consideration of ambient environmental/weather conditions as large temperature swings and high winds may invalidate this assumption as significant heat transfer with the environment may occur.In all cases, some heat losses are expected during the course of operation of the physical systems which will result in an under-sized design proportional to the actual heat losses, under these assumptions.
Additionally, heat losses from piping and heat exchangers are expected to be minimal.Similarly, these assumptions are valid for well-insulated systems.To account for any of these losses, the design engineer may apply higher-fidelity models that take into account operating temperatures, heat exchanger designs, and overall heat transfer coefficients for the system(s) of concern.However, doing so at this stage significantly complicates the model from an optimization perspective and negatively impacts computational tractability.Alternatively, high-fidelity modeling may be used to estimate average losses for a particular geographic region and IPH temperature requirements.These losses can then be accounted for as an overall thermal efficiency and used to adjust the IPH demand accordingly, similar to the factor 6 in the optical efficiency model for the CST array in Equation (1).Finally, to assess the performance of a CST system with respect to its ability to augment the use of conventional energy sources for IPH, we define the annual average solar fraction: In terms of the discrete values calculated from Equation ( 8) and the identities in Equation ( 9), the solar fraction can be approximated numerically using: which represents the fraction of total energy consumed by the process that is supplied by the solar system over the entire time interval of interest (e.g., a TMY).Lastly, although it was assumed that qp is not a function of the state or design decisions of the solar energy system, it need not be constant.For the purposes of studying the effects of transient load demand on the solar hybridization strategy, we will consider the model: which assumes i is in hours (h = 1) with i = 1 corresponding to 1:00 a.m. on 1 January, qp is the mean IPH demand, and σ p ∈ [0, 1] is the fractional deviation from the mean IPH demand.For σ p = 0, this models a constant IPH demand.For σ p > 0, this model accounts for a dynamic ramp-up and ramp-down of demand with a maximum demand at local noon and a minimum demand at local midnight.
To use the developed models effectively within a simulation-based optimization framework, we present the following concavity result for SF.Theorem 1.Let X a , X ts ⊂ R + be nonempty compact intervals.The solar fraction SF : X ts × X a → R is a concave function on its domain.
Proof.From Equations ( 7)-( 9) and the feasible system states, we can write for the qi a > 0 case (or else qi l = 0).The physical meaning of this case is that the solar system is producing more thermal power than can be used by the process or stored in the TES, and hence lost/rejected.We can see that h( qi s − qi p ) is an affine function of x a and is a convex function of x ts and x a .Therefore, is an affine function of x ts and x a .Since qi l is the pointwise maximum of an affine function and a constant, it is convex.From Equation (10), we have the summation of the terms: which is an affine function of x a subtracting a convex function of x ts and x a , and hence concave.Finally, since SF in Equation ( 10) is a summation of concave functions (the denominator is not a function of x ts or x a ), it is concave.

Smooth Approximation
Since the energy balance equations from the previous section involve min and max operators, differentiability of SF as defined in Equation ( 10) cannot be guaranteed.Since gradient-based optimization algorithms typically require all functions to be at least once-continuously differentiable, we must develop a smooth formulation.Consider the following operator: as a smooth approximation of the binary max operator: max{0, z}.The maximum error between the max operator and max is √ /2 and occurs at z = 0.For |z| >> , the error approaches zero.Hence, for small values, this approximation is not expected to introduce appreciable error into the model.A smooth approximation of SF, as defined in Equation (10), is given by: Again, the motivation of the model development is for use within simulation-based (global) optimization.The following result ensures that the smooth approximation SF s is concave and differentiable.
Theorem 2. Let X a , X ts ⊂ R + be nonempty compact intervals.Then, the smooth solar fraction SF s : X ts × X a → R is concave and continuously differentiable.
x ts for some i, and f : R → R : z → min {0, z}.It is clear that f is convex on R and g i is affine.Therefore, qi l = f • g i is convex for every i.Further, both f and g i are continuously differentiable on their domains.Therefore, qi l = f • g i is continuously differentiable on R 2 .Since SF is defined as the sum of the terms with a constant denominator (with respect to x ts and x a ), it is concave and continuously differentiable on R 2 .

Economic Model
An economic analysis is at the heart of all new-technology investment decision-making.Therefore, successfully designing a solar-IPH hybridization strategy requires a thorough economic analysis to ensure an optimal venture.Often, such an analysis consists of a comparison of conventional energy prices and the levelized cost of electricity (LCOE) for electrical power systems or the levelized cost of heat (LCOH) for thermal systems.For IPH applications, the LCOH analysis is often appropriate.However, the LCOH is defined at a high-level as:

LCOH ≡
sum of costs over lifetime sum of solar heat utilized , which cannot be guaranteed to be convex under appropriate assumptions.In contrast, we will develop the economic model of total lifecycle savings (LCS) which represents the total cost savings (which we wish to maximize) realized by hybridizing the IPH system with CSP.
The lifecycle savings is defined as the difference between the total lifetime cost of energy for a conventional system and that of the hybridized system, accounting for the time-value of money and energy-price inflation: where x = (x ts , x a ) is simply the vector-form of our decision variables, t li f e is the total lifetime of the project in years, r is the discount rate, C p,i is the cost of conventional fuel at year i, C k cap,i is the annualized cost of capital including debt service at year i for financing model k, and C om,i is the annual operating and maintenance (O&M) cost of the CST-TES system at year i.The solar hybridization design objective then is to maximize f k LCS , and hence, it is ideal if f k LCS is concave on its domain.The optimization problem formulation will be discussed in the next section.
The cost of conventional fuel at year i is given by: where C th,0 is the current specific thermal energy cost and r th is an inflation rate on the fuel price.The O&M cost is typically estimated under a fixed-pricing model [42] (i.e., cost per kWh of thermal energy produced).
The annualized capital expenditure and debt service of the CST-TES system at year i is given by [32]: which assumes monthly-compounded interest at an annual rate of r c , over the loan term t debt (in years).
The term C k cap,0 is the total capital expenditure (debt principal), which can be calculated using different pricing models k.The most common equipment pricing model has the form: which accounts for discounts in volume pricing (economies of scale).The factors C PTC,0 and C TES,0 are cost parameters which can be found in Table A2.The exponents of x a and x ts and the cost parameters are derived from various vendor quotes for commercially-available and custom-built equipment.Note, this model is concave in the variables (x ts , x a ).Since this term is subtracted within f disc LCS , concavity of f disc LCS cannot be guaranteed using this model, in general, which complicates solving the optimal design problem (i.e., global optimization is required).
Alternatively, a fixed-pricing (linear) model can be used: where C PTC is the specific installed cost of the solar array ($/m 2 aperture area) and C TES is the specific installed cost of thermal energy storage ($/kWh).Taking C PTC and C TES as a nominal specific cost for projects at a scale relevant to the design problem will yield results that are sufficiently accurate for initial design feasibility purposes but could lack the accuracy necessary for detailed design and economic analysis.Further, without a priori knowledge of the expected scale of a design, accurate cost estimates cannot be guaranteed with this model.
The following result details the construction of the convex envelope of the concave discounted equipment pricing model in Equation (14).The convex envelope is the tightest convex lower-bounding function that can be constructed for any nonconvex function on a domain of interest.The importance here is that it provides a rigorous lower bound on Equation ( 14) over its domain of interest while providing the best-possible convex approximation of Equation ( 14) over that domain.
. The convex envelope of Equation ( 14) on X is given by with f a (x a ) = C PTC,0 x 0.92 a and f b (x ts ) = C TES,0 x 0.91 ts .
Proof.Since C disc cap,0 is separable, we can write C disc cap,0 = f ts (x ts ) + f a (x a ).Further, its convex envelope can be constructed as the sum of the convex envelopes of the separable factors f ts and f a [43].The convex envelopes of these (univariate) factors are calculated simply as the secant lines between their endpoints: Proof.Since x L ts = x L a = 0, we have f ts (x L ts ) = f a (x L a ) = 0 and Equation ( 16) reduces to From the above results, we have the following relationships between the models: for f f ix LCS constructed via Corollary 1.

Optimization
The solar hybridization optimal design problem is formulated as a constrained nonlinear program (NLP) using the economic objective in Equation ( 13) developed in the previous section.The model is formulated as where g : R 2 → R n c are n c performance constraints that we may include, and X is simply a two-component interval vector with x L and x U as lower-and upper-bounds on the variable vector x = (x ts , x a ), respectively.In this work, we define only one constraint: where ξ represents some minimum fraction of the total IPH energy supply that must come from solar (e.g., a social responsibility constraint or a tax-incentive constraint).However, this general formulation can accommodate any relevant set of constraints.The parameters of the optimization model in Equation ( 19) can be found in Table A2.
Theorem 4. The constrained NLP in Equation (19) with the inequality constraint g(x) = ξ − SFs(x ts , x a ) ≤ 0 and k ∈ {cv, f ix} is a convex program.
Proof.Since SF s is concave (Thm.2) and C k cap,i is convex for k ∈ {cv, f ix} and for all i (Corollary 1), f k LCS , k ∈ {cv, f ix} is concave by design.Further, since SF s is concave, the function 19) is a concave maximization problem on a convex feasible set, it is a convex NLP.
Note that, from Equation ( 18), we have the following relationship between the optimal solution values: LCS constructed via Corollary 1.Therefore, we can use the fixed-pricing model in Equation ( 15) for initial high-level economic feasibility studies, which will overestimate the value of the more precise discount-pricing model in Equation (14).Further, given the result of Theorem 3, we can use Equation ( 19) with k = cv as a valid upper-bounding problem and solve the nonconvex program in Equation ( 19) with k = disc to global optimality by applying a spatial Branch-and-Bound algorithm (see Appendix C).Therefore, for greater-detail engineering design and economic feasibility studies, a global optimum of Equation ( 19) with k = disc can be readily obtained using the models presented herein.

Numerical Experiments and Results
The models presented above were implemented in the Julia programming language [44] using JuliaPro v.0.6.2.2 (Cambridge, MA, USA).Three disparate geographic locations across the US were considered as test cases for this model: Firebaugh, CA (California's Central Valley), Aurora, CO (Greater Denver Area), and Weston, MA (Greater Boston Area).The TMY hourly data for each location was obtained from the NREL NSRDB [41] (Physical Solar Model v.3).The monthly and annual average DNI values for each location can be found in Table A3.Note that the full TMY hourly dataset was used in this study.Two IPH demand profiles (i.e., constant and periodic) with equivalent total daily demand were considered for each geographic location, as were two fuel costs (i.e., commercial rate and industrial rate from [45]).For all studies, O&M cost was omitted (C om,i = 0) to simply illustrate the trade-off between conventional fossil energy consumption cost and solar capital costs.Additionally, for each study, the fixed-pricing model in Equation ( 15) calculated via Corollary 1 and the discount-pricing model in Equation ( 14) were used.
An optimal hybridization design study was conducted for each of the three geographic locations, for each case below.All problems were solved on a personal workstation running Windows 10 v1803 operating system with an Intel Core i7-5960X CPU and 32GB of RAM.All local optimization was performed using the interior-point algorithm IPOPT [46] from within the JuMP modeling platform [47] v0.18.0 utilizing forward-mode automatic differentiation for exact derivative information.Since f f ix LCS is concave on its domain, IPOPT is effective in solving Equation (19) with k = f ix (constructed using Corollary 1) to global optimality.Since f disc LCS is nonconvex on its domain, a global optimum cannot be guaranteed using IPOPT (or any local optimizer) alone.In this case, the problem is solved to global optimality using the the deterministic spatial Branch-and-Bound (B&B) algorithm in the EAGO software package [48,49] v0.1.0with a relative convergence tolerance of = 10 −2 and the greatest relative-width midpoint branching heuristic.The B&B algorithm for maximization is presented in Algorithm A1 of Appendix C with the concept of branching illustrated in Figure A1.The concave upper-bounding problem of Equation ( 19) with k = disc is simply given by k = cv.A valid upper bound is obtained at each iteration by solving the upper-bounding problem on the domain of interest using the IPOPT algorithm from JuMP.A valid lower bound is obtained by solving Equation (19) with k = disc to local optimality using IPOPT from JuMP.A global optimum is provided for each case for comparison against the optimal solution of the fixed-pricing model.

Constant IPH Demand, Commercial Fuel Rate
In this study, the IPH demand was modeled using Equation (11) with qp = 10 5 kW th and σ p = 0.
Hence, we have a constant IPH demand: qpeak p = qp = 10 5 kW th .The fuel cost was set at the commercial rate [45] in Table A2 and we let ξ = 0.
Table 1 contains the location information and corresponding optimal solar-IPH hybridization strategies (i.e., the optimal solutions of Equation ( 19)) for both the linear fixed-pricing model in Equation ( 15) and the nonconvex discount-pricing model in Equation (14).For each location, the fixed-pricing model determines an optimal design that has a lower TES capacity and less PTC aperture area (and therefore lower annual solar fraction) than the global optimal solution.The fixed-pricing model overestimates the total lifecycle savings by 6.5% for Weston, 1.7% for Aurora, and 1.6% for Firebaugh.For initial feasibility studies, this overestimation is sufficiently small as the optimal TES capacity and PTC aperture area are within roughly 5% of the global optimal designs for the discount-pricing model.
Table 1.The optimal solar-IPH hybridization strategies for the three geographic locations for constant IPH demand with qp = 10 5 kW th for the commercial fuel rate.The global optimal results for the concave capital cost model in Equation ( 14) are included for comparison against the linear capital cost model in Equation (15).The monthly thermal performance profiles for the global optimal designs are shown in Figure 3.Both Aurora and Weston experience less seasonal variance in DNI than Firebaugh (see Table A3) and so the thermal performance profiles appear flatter over the year than those for Firebaugh.Even though the annual average DNI for Aurora and Firebaugh are within about 8% of each other, they differ by as much as 36% in the summer months.Hence, the optimal designs for Firebaugh are significantly smaller and provide significantly more lifecycle savings than those for Aurora.Interestingly, the optimal TES capacity for Weston and Firebaugh are very similar, whereas the PTC aperture area for Weston sits between Firebaugh and Aurora.Clearly, the widely differing DNI profiles has a dramatic influence on the optimal designs.The f f ix LCS and SF s functions are plotted in Figure 4 for each geographic location for comparison.Each of the surfaces are qualitatively similar but there are clear differences between the geographic regions in terms of how sharply the economics change as you move away from the optimal design.The lifecycle savings for Weston are much more sensitive to the design decisions than Aurora and Firebaugh, which exhibit broader relatively-flat regions around the optimum.This explains why the optimal economics between the fixed-pricing and discount-pricing models are within less than 2% for Aurora and Firebaugh but almost 7% for Weston.It can be concluded from this analysis that designing solar-IPH hybridization strategies for lower-DNI regions may require the use of higher-accuracy cost models (e.g., discount-pricing models) with global optimization.

Periodic IPH Demand, Commercial Fuel Rate
In this study, the IPH demand is modeled using Equation (11) with qp = 10 5 kW th and σ p = 0.1.Hence, this case ramps up thermal power demand by 10% over the mean at the height of the workday ( qpeak p = 11 MW th ) and then ramps down to 10% below the mean by midnight to simulate reduced staffing and production output overnight.We consider the commercial rate for fuel [45] and we let ξ = 0.
Table 2 contains the location information and corresponding optimal solar-IPH hybridization strategies (i.e., the optimal solutions of Equation ( 19)) for both the linear fixed-pricing model in Equation ( 15) and the nonconvex discount-pricing model in Equation (14).As expected, the optimal designs have lower TES capacities than the constant IPH demand since the periodic demand is relatively synchronized with the sunrise and sunset.Further, the optimal PTC aperture areas are very similar to the constant IPH demand case since the total energy consumption between the two cases are the same.Note, the TES capacity x ts is scaled against the peak IPH demand qpeak p .Hence, to compare the TES capacity between the periodic IPH case and the constant IPH case, we need to scale them using qpeak p , for each case.Thus, for the Weston fixed-pricing model and constant IPH demand, we have a TES capacity of 11.24 h ×10 5 kW th =1.124 MWh th .For the Weston fixed-pricing model and periodic IPH demand, we have a TES capacity of 9.537 h ×1.1• 10 5 kW th =1.05 MWh th .Thus, the TES capacity is roughly 7% larger for the constant IPH demand over the periodic IPH demand.The periodic IPH demand is clearly more economical from a solar hybridization perspective.However, without accounting for the economics of the process itself, one cannot conclude that it should be operated in such a periodic fashion to achieve overall-better economics.In other words, it may be more costly to operate the process in a periodic fashion than what is saved with an optimal solar hybrid design.
Table 2.The optimal solar-IPH hybridization strategies for the three geographic locations for periodic IPH demand with qp = 10 5 kW and σ p = 0.1 for the commercial fuel rate.The global optimal results for the concave capital cost model in Equation ( 14) are included for comparison against the linear capital cost model in Equation (15).

Constant IPH Demand, Industrial Fuel Rate
We consider the IPH demand modeled using Equation ( 11) with qp = 10 5 kW th and σ p = 0 (i.e., constant demand profile).In this case, the industrial rate for fuel [45] is considered from Table A2 and we let ξ = 0.The objective function surfaces for both the fixed-pricing model (k = f ix) and discount-pricing model (k = disc ) are plotted in Figure 5. Table 3 contains the location information and corresponding optimal solar-IPH hybridization strategies (i.e., the optimal solutions of Equation ( 19)) for both the linear fixed-pricing model in Equation (15) and the nonconvex discount-pricing model in Equation (14).
Table 3.The optimal solar-IPH hybridization strategies for the three geographic locations for constant IPH demand with qp = 10 5 kW for the industrial fuel rate.The global optimal results for the concave capital cost model in Equation ( 14) are included for comparison against the linear capital cost model in Equation (15).It is clear from Table 3 and Figure 5 that the lower fuel costs dramatically reduce the lifecycle savings, and therefore jeopardize the economic viability of the hybridization project.Further, the reduced fuel costs ensure that the capital cost terms have a much greater impact on the lifecycle savings.As a result, significant nonconvexity is observed for the discount-pricing model.In fact, for the Aurora and Firebaugh cases, there exists a suboptimal local minimum at x = x L = (x L ts , x L a ) ≈ 0. IPOPT was applied to the nonconvex problem directly and only found the suboptimal local minima unless the domain was manually partitioned or bounded.Hence, without the application of the B&B algorithm (Algorithm A1), one may falsely conclude that solar-IPH hybridization is not economically viable.For the Weston, MA case, the fixed-pricing model concludes that the project is not economically viable, and hence, neither is the project with the discount-pricing model due to the upper-bounding relationship in Equation (18).Similar trends in the optimal solutions for the other locations are observed as for the previous cases.However, it is worth noting that for the Aurora and Firebaugh cases, the economics between the fixed-pricing and the discount-pricing models are significant.This is because significant nonconvexity is introduced for the lower fuel cost case.Hence, one can conclude that, for this low fuel cost case, the fixed-pricing model calculated via Corollary 1 over this large of a design space will not yield sufficiently accurate results.As a result, the discount-pricing model should be used with deterministic global optimization to ensure accurate results are obtained for feasibility studies.

Periodic IPH Demand, Industrial Fuel Rate
We consider the IPH demand modeled using Equation (11) with qp = 10 5 kW th and σ p = 0.1 (i.e., fluctuating demand profile).In this case, the industrial rate for fuel [45] is considered from Table A2 and we let ξ = 0.
Table 4 contains the location information and corresponding optimal solar-IPH hybridization strategies (i.e., the optimal solutions of Equation ( 19)) for both the linear fixed-pricing model in Equation ( 15) and the nonconvex discount-pricing model in Equation (14).The results are similar to the constant IPH demand case for the industrial fuel rate.Global optimization is required to locate the global optimal solution for the discount-pricing model case.Interestingly, the optimal solar design for the periodic IPH demand case is measurably larger in PTC aperture area over the constant IPH demand case for the same industrial fuel rate.This result differs from the commercial fuel rate cases which showed nearly identical aperture sizes between the designs.The results we see here is largely due to the very low annual solar fraction and the fact that the optimal designs call for the minimum TES capacity.Hence, the periodic IPH demand enables a greater solar fraction without adding TES capacity.However, as was identified for the constant IPH demand case, the fixed-pricing model dramatically overestimates the lifecycle savings due to the large degree of nonconvexity introduced by the discount-pricing model in the lower fuel cost case.As a result, it is still recommended for this case that the discount-pricing model is used with deterministic global optimization to ensure accurate results are obtained for feasibility studies.Table 4.The optimal solar-IPH hybridization strategies for the three geographic locations for periodic IPH demand with qp = 10 5 kW for the industrial fuel rate.The global optimal results for the concave capital cost model in Equation ( 14) are included for comparison against the linear capital cost model in Equation (15).

Computational Efficiency
The solution times for each of the numerical experiments are provided in Table 5.The computational cost trade-off of using the fixed-pricing model versus the discount-pricing model is clear.In almost all cases, using the discount-pricing model requires an order-of-magnitude longer to solve the problem than the fixed-pricing model due to nonconvexity and the need for the B&B algorithm.Despite this, the solution times for the discount-pricing model are still favorable considering the model's scale with O(10 4 ) state variables.As previously mentioned, the Weston location with the lower industrial fuel rate (for both IPH demand profiles) yielded the result that the hybridization project would not be economically viable.As can be seen in Table 5, these are the only cases where solving the discount-pricing model requires roughly the same order of computational effort as solving the fixed-pricing model.Additionally, due to the relationship in Equation ( 18), one can immediately conclude from the solution of the fixed-pricing model that the project is not economically viable.Therefore, for the purposes of investment decision-making, an additional "economic feasibility" termination criterion can be added to Algorithm A1, whereby the algorithm simply terminates if the upper bound is below some threshold of economic viability.For the purposes of this optimal design problem, the optimal solution translates to not implementing any CSP, and hence, a termination criterion could also be added that identifies this case.With these termination criteria, the B&B algorithm (Algorithm A1) would be identical in solution time to the fixed-pricing model since it would terminate at the root node with the solution of economic infeasibility after solving only the upper-bounding problem.For such cases, A2, the lower bounds on the decision variables were not exactly zero, albeit very close with respect to their scale, which helps avoid some numerical round-off issues observed.Since x L ≈ 0, the results still hold.

Conclusions
A dynamic concentrating solar power performance model is presented for use with design optimization and investment decision-making problems for hybridizing IPH systems.The model utilizes high-fidelity solar resource data from the NSRDB [41] and accounts for the aperture area of a parabolic trough solar concentrator and the storage capacity of a thermal energy storage device as design decision variables.The solar performance model was constructed such that it is concave and differentiable on the domain of interest.Since determining the performance of the solar system requires a dynamic simulation (with 8760 time steps using hourly DNI data), concavity of the model on the decision domain ensures that it is efficient for use with mathematical optimization.Additionally, an economic model based on the lifecycle savings was constructed that utilizes the solar performance model as well as two capital cost models: a linear fixed-pricing model and a nonconvex discount-pricing model.The convex envelopes of the nonconvex discount-pricing model were presented, which, when used within a spatial Branch-and-Bound algorithm, enable the global solution of the economic optimization problem with the nonconvex discount-pricing model.For the higher-cost commercial fuel rate, optimization of the fixed-pricing model and the discount-pricing model show good agreement.However, for the lower-cost industrial fuel rate, significant nonconvexity complicates the problem and therefore the greater-accuracy discount-pricing model (and global optimization) is recommended for investment decision-making.
The economic model presented herein is relatively simple to illustrate the trade-off between the fuel savings and the capital cost of the project.Greater levels of detail can easily be incorporated in this model, such as taxes, incentives, capital depreciation, carbon pricing, detailed O&M, etc.Similarly, the solar thermal performance model assumes the simplifying characteristics of negligible heat losses, ideal temperature delivery, and negligible CST array shadowing.These details can be readily incorporated in this model, if necessary.The developed models may serve as the foundation for more greater-complexity model-building for systems performance and economic analyses.Specifically, the models may easily incorporate various sources of uncertainty and investment and/or thermal performance thresholds to address robust design and operations under uncertainty and new technology investment risk.However, if enhancing complexity, the user must consider how convexity/concavity and computational tractability are affected.For example, simple tax models and carbon pricing models may be incorporated as increased O&M costs without significantly affecting computational cost or convexity/concavity results.However, high-fidelity heat transfer models complicate the calculation of the annual solar fraction SF s which substantially increase the integration cost and may introduce nonconvexity/nonconcavity. In this case, solving the optimization problem(s) will become extremely computationally expensive as convex and concave relaxations of SF s must be introduced, ensuring many more iterations of the B&B algorithm are required to identify an optimal solution.Similarly, accounting for various sources of uncertainty may have equivalent effects.As additional modeling complexity may add significant computational burden, the models presented herein represent a balanced approach between computational cost and modeling complexity for the purposes of optimization problem tractability and sufficient accuracy of the results for engineering design and investment decision-making.Definition A4 (Zenith Angle, φ [27]).The solar zenith angle φ (in degrees) is the angle between the sun's rays and the vertical plane: Definition A5 (Solar Elevation (Altitude) Angle, α [27]).The solar elevation or altitude angle α (in degrees) is the angle between the sun's rays and the horizontal plane: For certain periods of the year and for certain locations, this may not be the case (i.e., cos γ ≤ tan δ/ tan[Latitude • ]), and the azimuth angle becomes: Definition A7 (Surface Tilt Angle, β [27]).The surface tilt angle β of a horizontal surface with a N-S axis orientation and E-W tracking is given by: Definition A8 (Surface Incidence Angle, θ [25,27,50]).The surface incidence angle θ (in degrees) is the angle between the sun's rays and the normal of the horizontal surface.For fixed N-S axis with E-W tracking, the incidence angle can be calculated as: cos θ = cos 2 φ + cos 2 δ sin 2 γ.

Figure 1 .
Figure 1.(left) An illustration of a parabolic trough solar concentrator (PTC) with the defined array length L m , aperture length L a and the aperture plane with incident solar radiation.(right) The large-aperture PTC (Skyfuel, Lakewood, CO, USA) powering a solar thermal desalination pilot by Stuber et al. [38].

Figure 2 .
Figure 2. The block-flow diagram representation of the hybrid solar industrial process heat system being modeled.

Figure 3 .
Figure 3.The thermal performance of the global optimal hybridization strategies for constant IPH demand and the commercial initial specific fuel cost.

Figure 4 .
Figure 4.The objective function surfaces f f ix LCS (left column) and the solar fraction surfaces SF s (right column) for constant IPH demand for commercial fuel rate show similar qualitative trends between each region (rows).

Figure 5 .
Figure 5.The objective function surfaces f f ix LCS (left column) and f disc LCS (right column) for constant IPH demand for industrial fuel rates are plotted for each geographic location.The discount-pricing model exhibits nonconvexity in all cases and suboptimal local minima for Aurora and Firebaugh.
with B in degrees and ST in minutes and the integer N as the day number.To relate the solar time to the local time, a time correction, TC (in minutes), is required:TC ≡ ST ± 4[15(Time Zone (h)) − Longitude • ]where (+) is chosen if the location is east of the prime meridian and (−) if it is to the west.Definition A2 (Declination Angle, δ[27]).The solar declination angle δ is the angular distance of the sun's rays north (or south) of the equator: δ ≡ 23.45 • sin B. Definition A3 (Hour Angle, γ [27]).The hour angle γ is the angle between the meridian of an observer and the sun.It is negative before local solar noon and positive after local solar noon.At local solar noon, γ = 0 • .The hour angle (in degrees) is defined as γ ≡ 15(Local Time Hour + TC/60 − 12).

Table 5 .
The solution times (in seconds) for each of the numerical experiments.