Optimal Switchable Load Sizing and Scheduling for Standalone Renewable Energy Systems

TL;DR

This paper presents analytical solutions and algorithms to optimize load sizing and switching schedules in off-grid renewable systems, maximizing solar energy utilization without costly storage, demonstrated with real solar data.

Contribution

It introduces novel mixed-integer linear programming and constrained least squares algorithms for optimal load sizing and switching in off-grid systems, enhancing solar energy utilization.

Findings

73% solar energy utilization with two loads

98% utilization with six loads

Algorithms outperform baseline scheduling methods

Abstract

The variability of solar energy in off-grid systems dictates the sizing of energy storage systems along with the sizing and scheduling of loads present in the off-grid system. Unfortunately, energy storage may be costly, while frequent switching of loads in the absence of an energy storage system causes wear and tear and should be avoided. Yet, the amount of solar energy utilized should be maximized and the problem of finding the optimal static load size of a finite number of discrete electric loads on the basis of a load response optimization is considered in this paper. The objective of the optimization is to maximize solar energy utilization without the need for costly energy storage systems in an off-grid system. Conceptual and real data for solar photovoltaic power production is provided the input to the off-grid system. Given the number of units, the following analytical solutions and computational algorithms are proposed to compute the optimal load size of each unit: mixed-integer linear programming and constrained least squares. Based on the available solar power profile, the algorithms select the optimal on/off switch times and maximize solar energy utilization by computing the optimal static load sizes. The effectiveness of the algorithms is compared using one year of solar power data from San Diego, California and Thuwal, Saudi Arabia. It is shown that the annual system solar energy utilization is optimized to 73% when using two loads and can be boosted up to 98% using a six load configuration.

Tables

Table 1. Table 1: ECLS results for unit sizes and solar utilization for n 𝑛 n between 2 and 6.

$n$	$x_1$	$x_2$	$x_3$	$x_4$	$x_5$	$x_6$	$\sum \left(x\right)$	$SU$
2	0.4670	0.1450					0.6120	0.7071
3	0.5068	0.2025	0.0857				0.7950	0.8277
4	0.5111	0.2269	0.1175	0.0585			0.9140	0.9117
5	0.4769	0.2138	0.1117	0.0565	0.0290		0.8880	0.9571
6	0.4476	0.2015	0.1060	0.0544	0.0287	0.0158	0.8540	0.9790

Table 2. Table 2: Same as Table 1 , but for ICLS method.

$n$	$x_1$	$x_2$	$x_3$	$x_4$	$x_5$	$x_6$	$\sum \left(x\right)$	$SU$
2	0.4078	0.1994					0.6053	0.7274
3	0.4210	0.2076	0.1028				0.7314	0.8601
4	0.3957	0.1954	0.0989	0.0467			0.7367	0.9273
5	0.4180	0.2063	0.1034	0.0508	0.0228		0.8013	0.9614
6	0.3913	0.1935	0.0973	0.0473	0.0233	0.0115	0.7642	0.9796

Table 3. Table 3: Same as Table 1 , but for Mixed-Integer Linear Programming.

$n$	$x_1$	$x_2$	$x_3$	$x_4$	$x_5$	$x_6$	$\sum \left(x\right)$	$SU$
2	0.3933	0.1899					0.5832	0.7274
3	0.4126	0.1899	0.1008				0.7033	0.8586
4	0.4074	0.2100	0.0914	0.0438			0.7526	0.9243
5	0.3811	0.2150	0.1082	0.0515	0.0240		0.7798	0.9597
6	0.4505	0.2049	0.0837	0.0598	0.0327	0.0162	0.8478	0.9765

Table 4. Table 4: Comparison of the advantages and disadvantages of the optimization algorithms.

	ICLS	ECLS	MILP	Analytic
Advantage	-Fast -Considering all sorted data -Not equidistant	-Fast -Unique minimum for fixed constraint	-Convex (Global optimality)∗ -Any data profile for limited variables - More constraints could be added	-Fast -Exact (no approx.) -Guarantees global optimum
Disadvantage	-No guaranteed global optimality -Sensitive to initial conditions	-Down-sampling required -Equidistant -Requires line-search for constraint	-Computationally expensive -Requires down-sampling	-Only for symmetric, i.e., clear days -Can not be used for planning algorithms
∗ Note: convexity for MILP approach is assumed in case of convergence

Equations

$$S_{max}=\underset{k}{\max }S(t_k),$$

$$E(t_k)=S(t_k)-\sum _{i=1}^n{u}_k^i{x}^i$$

$$\begin{array}[]{rcll}u_{k}&=&\left[\begin{array}[]{cccc}u_{k}^{1}&u_{k}^{2}&...&u_{k}^{n}\end{array}\right],{}&{u_{k}^{i}}\in[0,1]\\ x&=&\left[\begin{array}[]{cccc}x^{1}&x^{2}&...&x^{n}\end{array}\right]^{T},{}&x^{i}>0\end{array}$$

$$E(t_k)=S(t_k)-u_{k}x,$$

$$\text{mbox}arg\underset{u_k,x}{\min }\sum _{k=1}^{N}E(t_k{)}^2,E(t_k)=S(t_k)-u_{k}x$$

$$x^i\ge x^{i-1}\ge 0$$

$$\sum _{i=1}^n{x}^i=|x{|}_1=C{S}_{max}$$

$$A^{\prime }(t)=2S(t)+2t{S}^{\prime }(t),\text{or}{A}^{\prime }(y)=2S^{-1}(y)+2y(S^{-1}(y){)}^{\prime }$$

$$S^{-1}(y)=-\bar{y}(S^{-1}(y){)}^{\prime }$$

$$GHI=951.39cos(t{)}^{1.15}$$

$$S(t)=y=acos(bt)+a_{1}sin(bt)=asin(bt+c)$$

$$t=S^{-1}(y)=\alpha \arcsin (\beta y)+\gamma$$

$$t^{\prime }=(S^{-1}(y){)}^{\prime }=\frac{\alpha \beta }{\sqrt{1+{\beta }^2{y}^2}}$$

$$S^{-1}(y)+y(S^{-1}(y){)}^{\prime }=0\to \frac{y}{\alpha \beta \sqrt{1-y^2/{\alpha }^2}}-\frac{\gamma -arcsin(y/\alpha )}{\beta }=0$$

$$\begin{array}[]{rcl}S(\bar{t})=\bar{y}&=&0.6489\\ \bar{t}&=&\pm 123\end{array}$$

$$\text{SU}=\text{Solar utilization}=\frac{\text{Total Energy Captured by Units}}{\text{Total Solar Energy}}=\frac{2\bar{t}\bar{y}}{{\sum }_{-t_{max}}^{t_{max}}S(t)dt}=56.03\%.$$

$$A(y_1,y_2,y_3)=\sum _{i=1}^{2^n-1=3}{A}_i=2S^{-1}(y_2)y_2+2(S^{-1}(y_1)-S^{-1}(y_2))y_1+2S^{-1}(y_3)(y_3-y_2)$$

$$\bar{A}(y_1,y_2)=2y_1(S^{-1}(y_1)-S^{-1}(y_2)+S^{-1}(y_1+y_2))+2y_2{S}^{-1}(y_2).$$

$$\displaystyle\begin{array}[]{c}\underset{{y}_{1},{y}_{2},{y}_{3}}{\text{max }}{A}({y}_{1},{y}_{2},{y}_{3})\\ \begin{array}[]{ll}\text{subject to}&0\leq{y}_{1}\leq{y}_{2}\leq{y}_{3}<y_{max}\\ &{y}_{3}={y}_{1}+{y}_{2},\\ \end{array}\end{array}$$

$$\displaystyle\begin{array}[]{c}\underset{{y}_{1},{y}_{2}}{\text{max }}{\bar{A}}({y}_{1},{y}_{2})\\ \begin{array}[]{ll}\text{subject to}&0\leq{y}_{1}\leq{y}_{2}<y_{max}\\ \end{array}\end{array}$$

$$▽J=\left[\begin{array}{c}\frac{\partial A}{\partial y_1}\\ \\ \frac{\partial A}{\partial y_2}\end{array}\right]=\left[\begin{array}{c}{S}^{-1}(y_1)-S^{-1}(y_2)+S^{-1}(y_1+y_2)+y_1(S^{-1}(y_1){)}^{\prime }\\ \\ (S^{-1}(y_2){)}^{\prime }(y_1+y_2)+S^{-1}(y_2)\end{array}\right]$$

$$\text{SU}=\frac{A({\bar{y}}_1,{\bar{y}}_2,{\bar{y}}_3)}{{\sum }_{-t_{max}}^{t_{max}}S(t)dt}=79.49\%$$

$$A(y_1,\cdots ,y_{2^{n-1}})=\sum _{i=1}^{2^n}2[S^{-1}(y_i)-S^{-1}(y_{i+1})]y_i,$$

$$y_i=bin(i{)}^T\times \left[\begin{array}{c}{y}_1\\ y_2\\ ⋮\\ y_n\end{array}\right],i=\{1,2,\cdots ,2^n-1\},$$

$$S(t_{\bar{k}+1})\ge S(t_{\bar{k}}),\bar{k}=1,...,T-1$$

$$S(t_{\bar{k}+1})\ge S(t_{\bar{k}})\Rightarrow u_{\bar{k}+1}x\ge u_{\bar{k}}x,\bar{k}=1,...,T-1$$

$$f(u_{\bar{k}})=\sum _{i=1}^n{f}^i{u}_{\bar{k}}^i,f^i>0$$

$$f(u_{\bar{k}+1})\ge f(u_{\bar{k}}),\bar{k}=1,...,T-1$$

$$f(u_{\bar{k}})=d_{\bar{k}}=\sum _{i=1}^n{u}_{\bar{k}}^i\cdot 2^{n-i},$$

$$\begin{array}[]{c}\displaystyle\hat{x}=\mbox{arg}\min_{x}\sum_{\bar{k}=1}^{T}E(t_{\bar{k}})^{2},~{}~{}E(t_{\bar{k}})=S(t_{\bar{k}})-u_{\bar{k}}x\\ f(u_{\bar{k}+1})\geq f(u_{\bar{k}}),~{}\bar{k}=1,...,T-1\end{array}$$

$$\begin{array}[]{rcl}S&=&\left[\begin{array}[]{cccc}S(t_{1})&S(t_{2})&\cdots&S(t_{T})\end{array}\right]^{T}\in R^{T\times 1},\\ E&=&\left[\begin{array}[]{cccc}E(t_{1})&E(t_{2})&\cdots&E(t_{T})\end{array}\right]^{T}\in R^{T\times 1},\\ U&=&{\left[\begin{array}[]{cccc}u^{1}_{1}&u^{2}_{1}&\cdots&u^{n}_{1}\\ u^{1}_{2}&u^{2}_{2}&\cdots&u^{n}_{2}\\ \vdots&&\ddots&\vdots\\ u^{1}_{2^{n}-1}&u^{2}_{2^{n}-1}&\cdots&u^{n}_{2^{n}-1}\end{array}\right]\otimes\mathbf{1}_{1\times L}}\in R^{T\times n}\\ x&=&\left[\begin{array}[]{cccc}x^{1}&x^{2}&\cdots&x^{n}\end{array}\right]^{T}\in R^{n\times 1}\end{array}$$

Full text

††journal: Solar Energy

Nomenclature

[TABLE]

1 Introduction

Increasing global energy demand and human population growth have triggered a need for standalone renewable applications. Recent estimates show that 1.4 billion people do not have access to energy services and one billion are suffering from unreliable electricity services [1]. Standalone application of clean energy, (E.g., fresh water pumping), has become more critical for humanity [1, 2]. Often such systems are powered by solar photovoltaic (PV) due to ubiquitous high solar resource availability and scalability. However, solar production exhibits high variability over a broad range of time scales [3]. Power variability is the main obstacle facing solar energy in standalone or islanded mode applications. High penetrations of solar power sources create large power swings which influence electric power quality [4], and can cause loss of load or generation curtailment [5]. Variability of solar PV generation is a result of seasonal and diurnal changes in the sunpath as well as short-lived cloud cover. Solar variability limits the operation of off-grid loads at maximum capacity [6, 7, 8].

Optimal load switching can be applied to microgrids with any hybrid forms of renewable energy resources such as solar and wind [9, 10, 11, 12] to capture as much renewable energy as possible. Although partial or modulated load operation is conducive to the problem, there are numerous types of load units which can only be switched on or off, such as non-dimmable lighting, standard electric motors, and Magnetic Resonance Imaging (MRI) machines at hospitals and load aggregation such as demand side management [13, 14]. Dispatching such binary load units, which are referred as switchable loads hereafter, to follow available renewable energy resources have been discussed in the literature for different microgrid applications such as water desalination [15], pumping systems [16], irrigation systems [17], and cooking appliances[18, 19, 20].

Different optimization techniques are used for planning and design of such systems. For instance, mixed-integer linear programming (MILP) has been used in many fields, such as unit commitment of power production [21] and power transmission network expansion [22, 23], as well as scheduling problem of the generation units in off-grid in order to maximize supply performance of the system [24]. Nonlinear approaches have also been applied to load scheduling [25]. For example, neural networks and genetic algorithms have been applied to size stand-alone PV [26, 27, 28]. The on/off control optimization problem is similar to the unit commitment problem in power systems and bio-fuel[29, 30]. However, limited studies have been conducted on the optimal load sizing in a standalone (islanded) grid application with switchable loads. Most of other research has been in the demand/supply side while very few looked into the unit/load sizing for many reasons, such as, the load is assumed to be fixed and has to meet by any supply way [31] or the accessibility of designing load for certain application is harder and not easy process. This work focuses on optimal load sizing for standalone applications in rural areas or off-grid sites. While the present paper assumes an off-grid system, similar challenges exist for a power system with a weak grid connection, i.e. with a line carrying capacity that could only balance variability that is a small fraction of local solar generation or load capacity.

Energy storage systems (ESS) have been applied to solve the variability challenges [32, 33].An alternative or complementary approach is optimal sizing and scheduling of load units which follow power generation variability to maximize solar energy utilization and load uptime. Clearly, the solar energy utilization could be improved with an ESS, but an ESS that eliminates solar variability would need to be large enough to store several days’ worth of solar power which is uneconomical at present. Smaller ESS would experience significant cycling and deep discharge events if not properly maintained, increasing maintenance costs and requiring replacement much before the end-of-life of a PV system. Our objective is to improve solar utilization without an ESS and use load demand response only and show that high efficiencies can still be obtained. In practice, a combination of a small ESS with high cycle life such as an ultracapacitor ESS and the proposed load sizing and scheduling system would probably be the best solution. The ESS would absorb solar variability at time scales of seconds to minutes while the loads would balance variability at longer time scales. This approach would allow limiting ESS energy capacity making it more economical. While practical challenges of implementing such a system are significant, e.g. in maintaining system stability during switching, this paper focuses on the critical algorithmic work that permits such a system to operate efficiently and economically.

In a properly planned system the solar system would be optimally sized to power the load required for the intended application. This paper does not consider this scenario. Often in practice the conditions are not as plannable. Load growth will occur and a solar power system may be initially oversized to accommodate such growth. Sizing the solar system may also be limited by land ownership and topographic constraints. The solutions proposed in this paper apply in such a context where solar capacity is fixed and loads are sized to optimize solar energy utilization.

This paper proposes an optimization model to capture the maximum amount of variable solar generation, which sizes and schedules a finite number of loads to track available solar PV power. The objective is to maximize solar utilization, given the projected power generation of the renewable energy resources. Here, solar utilization is defined as percentage of energy captured by the units over total solar energy produced. This is akin to terms such as solar utilization factor [34] and loss of power supply (LPS) [35] which are commonly used in the literature. The loads are assumed to switch between a binary ”on” or ”off” statuses, where both the switching times and the size of the static power demand (static load size) determines the ability to track available solar power.

The main contribution of this paper is to develop both analytical solutions and computational approaches based on Equality Constrained Least Squares (ECLS), Inequality Constrained Least Squares (ICLS), and Mixed-Integer Linear Programming (MILP) in order to solve the optimization problem. The rest of the paper is organized as follows. The mathematical formulation is given in Section 2 along with an analytical example and the motivation for a computational procedure for optimal load size selection. One year of solar resource data for San Diego is analyzed and discussed in Section3. Section 4 presents different computational procedures for optimal load size selection based on a bi-linear optimization problem involving a mix of binary and real numbers. The simulation results are presented and discussed in Section  5. Finally, Section 7 concludes the paper.

2 Problem Formulation

The sizing and scheduling problem computes the distribution of the optimal load size of a finite number load units, given an available (solar) power profile. For the optimal load size selection, the loads are assumed to operate in a binary manner, off or on, and therefore only the static load size is optimized.

2.1 Static Load Response Optimization Problem

To formalize the notation for the optimization approach presented in this paper, we assume knowledge of the (solar) power delivery $S(t)$ sampled at regular time intervals $t=t_{k}=k\Delta_{T}$, where $1/\Delta_{T}$ is a fixed sampling frequency and $k$ is the sample index. In this way, we have a data set of $T$ points on the solar power production $S(t_{k}),~{}k=1,...,T$. Typically, $S(t_{k})$ is close to a daily periodic function and $S(t_{k})\geq 0$ over a daily time interval $t_{k}\in[t_{b},t_{e}]$, where $b$ is the beginning and $e$ is the ending of the day, with a maximum value

[TABLE]

that is typically equal to the AC rating of the solar power system. The data $S(t_{k})$ and $S_{max}$ may be available from historical solar power measurements or from solar radiation measurements in conjunction with a model of the solar power conversion efficiency.

In the static load response optimization we consider $n$ loads, where each load $i=1,...,n$ is characterized only by a static power value $x^{i}$ that can be either turned on or off. Given the number $n$ of loads, the objective of the static load response optimization is to find the optimal distribution of static load values $x^{i},~{}i=1,2,\ldots,n$ so that the time sampled solar power delivery $S(t_{k})$ can be approximated as closely as possible to maximize the energy captured.

As indicated in Figure 1, the power mismatch $E(t_{k})$ at any time $t_{k}$ can be characterized by

[TABLE]

where $u_{k}^{i}\in[0,1]$ are $n$ binary numbers reflecting the on/off switch state of the individual loads $i=1,\ldots,n$ with their (to be determined) static load size $x^{i}$. Defining the vectors

[TABLE]

the static power mismatch $E(t_{k})$ at a particular time $t_{k}$ can be written with an inner product

[TABLE]

of the time dependent binary switch state vector $u_{k}$ and the static load size distribution vector $x$. Static load response optimization can now be written as

[TABLE]

where the variables $u_{k}$ and $x$ are given in (1). The optimization in (2) is a least squares optimization in which both the time dependent binary switch state vector $u_{k}$ and the static load size distribution vector $x$ must be determined on the basis of the $T$ data points on the solar production $S(t_{k}),~{}k=1,\ldots,T$.

Clearly, the least squares optimization in (2) is non-standard for several reasons. First of all, the error $E(t_{k})$ is bi-linear due to the product of the optimization variable $u_{k}$ and $x$. Furthermore, the optimization variable $u_{k}$ is a binary vector, whereas the elements $x^{i}$ of the static load size distribution vector $x$ must likely satisfy (linear) constraints

[TABLE]

to ensure a unique load distribution solution with real valued positive loads. An additional linear constraint

[TABLE]

where $C$ can be chosen in the range of ${0.5<C<1}$ ensures that the sum of the load distribution is bounded to avoid oversizing of the loads in trying to match the anticipated maximum power production $S_{max}$. Finally, the number of loads ($n$) also needs to be determined. It is clear that a larger value $n$ of will enable smaller power mismatch errors $E(t_{k})$ but would likely increase the investment cost as the per kWh cost of a load unit typically decreases with the size of the unit.

The process can be described in a flow chart as Figure 2 illustrates, starting with annual solar irradiance data for at least one year to capture the seasonal changes. Then followed by a power model to compute solar power and then sort solar power annual data followed by inputting number of desired units then the optimal unit sizing will be computed then followed by a daily unit scheduling.

3 Solar Resource Data and its Distribution

To capture the interannual variability of solar irradiance, solar resource data should be collected for several years, such as in the production of a typical meteorological year (TMY). In this paper, only one year of solar power generation is used to demonstrate the model application, but most large solar system developers rely on multidecadal modeled power production based on site adoption of long-term satellite records with short-term local measurements for their financial calculations [36]. Such long-term data would be preferable in practice although interannual variability of solar energy generation is small. For example, [37] specifies that the interannual variability of GHI for 7-10 years of measurement at Potsdam, Germany and Eugene, USA is about 5%.

Here, only one year of data was available which allows characterizing most of the important seasonal and diurnal variability.

3.1 Case Study 1: San Diego, USA

Multidecadal PV power projections are typically based on the measured and modeled GHI which is transposed to the direct, diffuse, and reflected radiation at the plane-of-array[38, 39, 40] and input into a PV performance model. We bypass the complexity in the PV power modeling by using a solar power generation dataset available at two sites. AC power production and GHI data were collected from a (91.6 kWDC and 100 kWAC) fixed tilt (non-tracking) polycrystalline PV system (Figure 3). The system was installed at 10 ∘ tilt and facing south, at the UC San Diego campus at $32^{\circ}53^{\prime}01.4"$N $117^{\circ}14^{\prime}22.6"$W. Data for one year, i.e., May 2011 through April 2012, were used and averaged over 15 minutes. The raw data is available at 1 s resolution and the algorithm can be applied to data at any temporal resolution. Nevertheless, switching loads over such short timescales is generally impractical and we assume instead that a small energy storage system modulates high frequency solar variability to create a supply that is stable over 15 min intervals. The solar power data for these specific sites are not symmetric over a day, and overcast conditions occur more frequently in the mornings.

The impacts of solar PV generation variability on the utilization of different combinations of discrete loads is illustrated through 2-dimensional (2D) histograms in Figure 4. In each subplot, similar to Figure 3, one year of solar power data is superimposed over one day and different numbers of units $n$ are used to track the solar power generation. The quantization idea of the solar power data $S(t_{k})$ for different numbers of units is illustrated. $N=$ $2^{n}-1$ is the number of discrete combinations of units for each case described in (1), which implies that with more units more discrete load levels can be served resulting in a better match with the solar generation. Qualitatively, the best combination of load sizes is expected to be the one that is able to best track the solar generation levels (or minimize power mismatch as in (2)) that are (i) large and (ii) occur frequently.

3.2 Case Study 2: Thuwal, Saudi Arabia

Another case was selected to prove the robustness of the algorithms to different input data. Thuwal, north of Jeddah, is a city located in the west coast of Saudi Arabia as shown in Figure 2 of [41]. Thuwal solar meteorology is predominantly clear and at a lower latitude and is therefore quite different from San Diego with days such as clear day (D1) in Figure 10 being more common. Data from a monocrystalline Silicon solar PV power plant at ($22^{\circ}18^{\prime}28.5"$N $39^{\circ}06^{\prime}17.1"$E) with tilt 20*∘* and azimuth of 133*∘* and 145*∘* (split in two different arrays) was collected by King Abdullah University of Science and Technology (KAUST).

4 Optimization Techniques

This section is initiated by an analytical solution for a small number $n\leq 2$ of loads is presented when the available power follows a symmetric daily profile. The analytical approach is followed by the motivation for a computational procedure to compute optimal load size distribution for a larger number $n\geq 2$ of loads when the available daily power profile is non-symmetric. After that three different optimization techniques are presented, discussed and compared.

4.1 Analytical and Motivating Example

Consider a (symmetric) time dependent power function $y=S(t)$ which must be followed by the rectangular power demands created by a simple on/off switching of a static load over a specified time period. System efficiency is optimized by finding the largest rectangular window (representing energy demand by the switchable load units) to be drawn under the power function $S(t)$. For a single load subjected to a symmetric $S(t)$, this problem reduces to selecting an optimal on/off switch time $\bar{t}$ to define the width of the rectangle $2\bar{t}$ and height $\bar{y}=S(\bar{t})$, as indicated in Figure 5.

The function $A(t)$ that parametrizes the area of the rectangle in Figure 5 is equal to $A(t)=2tS(t)$ and can be written in terms of $y$ as $A(y)=2S^{-1}(y)y$. The derivative or Jacobian of this function is given as

[TABLE]

By setting $A^{\prime}(y)=0$, the optimal value for the load size $\bar{y}=f(\bar{t})$ and the resulting switch time $\bar{t}$ can be found. In this case, it is clear that $A^{\prime}(y)=0$ leads to the analytical expression

[TABLE]

where solving it for $y$ gives the optimal load size $\bar{y}$.

To make this approach applicable to a symmetric solar PV power curve, a clear sky power solar function $S(y)$ can be modeled as a trigonometric function or a parabola. Global horizontal irradiance (GHI) models have been extensively studied in the literature in [42, 38, 39, 40] take for example equation 20 [42] as bellow

[TABLE]

In this paper, a simplified fitting curve is presented using trigonometric and parabola functions as below (see Appendix for accuracy discussion),

[TABLE]

where $a=0.9903,~{}a_{1}=-0.001192$, $~{}b=0.006952$, and $c=1.572$ based on data obtained from a clear solar day, where axes interception are $(0,\pm y_{max})=(0,0.9903)$ and $(\pm t_{max},0)=(\pm 226,0)$. The inverse function $S^{-1}$ is thus as follows,

[TABLE]

where $\alpha=143.8,~{}\beta=1.01,\text{ and }\gamma=-226.1$. More details on the computation of the numerical values of $a$, $b$ can be found in Appendix A. The first derivative of (7) can be expressed as

[TABLE]

Solving the analytic expression (5) numerically for an optimal integer value $\bar{t}$, with $S^{-1}(y)$ given in (7) and $(S^{-1}(t))^{\prime}$ given in (8) with a trigonometric function to compute the optimal values of $\bar{y}$ analytically, with the following numerical value of $A^{\prime}(y)=0$

[TABLE]

then leads to

[TABLE]

for the optimal switch time $\bar{t}$ and normalized load size $\bar{y}$.

The resulting optimal switch time and load size lead to a maximum (rectangular) area of $2\bar{t}\bar{y}=2\cdot 0.6489\cdot 123=159.6294$. With the known (symmetric) solar power curve $S(t)$, we can also compute $\sum_{-x_{max}}^{x_{max}}S(t)$ over the integer values $t\in[-226,226]$ to be 284.8962 and obtain

[TABLE]

This means that for a single load, the optimal rectangular area achieved under the given symmetric power curve $S(t)$ captures 56% of the total (solar) energy. Similar results can be obtained by a straightforward line search algorithm for the scalar value of $t$ as provided in Appendix B.

The approach of finding optimal on/off switching time and load size for a symmetric power curve can also be extended to the case of multiple loads. However, the solution is only analytically tractable for $n=2$ loads where the Jacobian becomes a two dimensional vector or a three dimensional vector with an additional equality constraint as indicated below.

Figure 6 illustrates the optimal load size distribution and switching to maximize solar utilization for $n=2$ loads. Under this scenario there exists 3 optimal switch times ${t}_{1}$, ${t}_{2}$ and ${t}_{3}$ for two (optimal) load sizes ${y}_{1}$, ${y}_{2}$, where ${y}_{3}={y}_{1}+{y}_{2}$ is used to indicate when both loads are on. The dark gray shaded area is due to switch time ${t}_{2}$, where the area is $A_{1}=2t_{2}y_{2}$. The area of the light gray shaded areas are functions of $t_{1}$, $t_{2}$ and $t_{3}$, where the area of the light gray rectangular areas is the sum of $A_{2}=2(t_{1}-t_{2})y_{1}$ and $A_{3}=2t_{3}(y_{3}-y_{2})$. The objective is to maximize the sum of the shaded areas

[TABLE]

to solve for the optimal switch times and load sizes as shown in Figure 6. The area ${A}({y}_{1},{y}_{2},{y}_{3})$ can be reduced to a function of only two variables ${\bar{A}}({y}_{1},{y}_{2})$ by substituting ${y}_{3}={y}_{1}+{y}_{2}$ to obtain

[TABLE]

The maximum solar utilization can now be expressed as an optimization problem:

[TABLE]

where $S(y)$ is given in (7), and $y_{max}$ is the $y$ intersection of $S(y)$. We consider only the positive values $y>0$ of the symmetric trigonometric approximation $S^{-1}(y)$ as indicated in Figure 6.

Since the optimization problem has affine equality constraints, its solution set is convex. Moreover, the objective function is a concave function. Thus, it has a single maximum or minimum $1/{A}({y}_{1},{y}_{2},{y}_{3})$ as indicated in Figure 7. As a result, solving the optimization in (11i) by an iterative gradient method will lead to the global maximum solution. The following Jacobian matrix, which is derived in Appendix C,

[TABLE]

can be used in an iterative gradient based method, leading to the optimal load size solutions $[\bar{y}_{1},\bar{y}_{2},\bar{y}_{3}]=[\bar{y}_{1},\bar{y}_{2},\bar{y}_{1}+\bar{y}_{2}]=[0.2727,~{}0.5758,~{}0.8485]$. The resulting solar utilization for 2 loads is then characterized by

[TABLE]

indicating a significant improvement over the single load solar utilization of $56.0307\%$ over the same symmetric solar power curve.

The analytic approach indicates that maximizing solar utilization is equivalent to finding the largest sum of rectangle windows that can be drawn under the power function $S(t)$. Extending this concept to $n$ number of loads where $n>2$ would entail

[TABLE]

where $bin(i)$ is a reversed vertical vector format representing the binary value of $i$. The number of variables of the new area function $\bar{A}({y}_{1},\cdots,{y}_{n})$ in (12a) reduces to $n$ variables instead of $2^{n}-1$ by substituting (12b).

In summary the optimization problem appears to be convex over our solution set. However, relying on $S^{-1}(y)$ for $y>0$ is not possible, as $S^{-1}(y)$ is not guaranteed to exist. Furthermore, a power curve $S(t)$ may not be symmetric, especially for PV systems operating on non-clear day conditions. To overcome these obstacles, this paper proposes optimization approaches that exploit the convexity of the optimization problem that selects the optimal load size distribution. Although the analytic approach is only viable for a small number of loads under symmetric power curves, we will use the optimally computed solar utilization as a benchmark for the solar utilization obtained from the optimization approaches presented in the following.

4.2 Equality Constrained Least Squares Optimization

Although the optimization in 2 to minimize the Least Squares of the static power mismatch $E(t_{k})$ is bi-linear, it is clear that the entries of the time dependent binary switch state vector $u_{k}$ is given by a limited number of binary combinations. The number of binary combinations depends on the choice of $n$ and the number of data points $T$. Once the time dependent binary switch state vector $u_{k}$ is fixed, the optimization in (2) reduces to a standard Least Squares (LS) problem to find the optimal value of the static load size distribution vector $x$.

To ensure a unique solution for the static load size distribution vector $x$, the constraint (3) can be included in the LS optimization implicitly by simply ordering the $T$ data points of the (solar) power data $S(t_{k})$. The ordering uses the fact that both $S(t_{k})\geq 0$, $x^{i}>0$ and the fact that a larger value of $S(t_{k})$ would require the switching of a larger sum of loads ${u_{k}}x$. To set up the solution to the optimization to (2), the solar data $S(t_{k})$ that may be periodic due to daily patterns and irregular due to weather patterns, is sorted such that

[TABLE]

The ordering in (13) ensures that $S(t_{\bar{k}})$ is monotonically non-decreasing function represented as the unshaded curve in Figure 8 In addition, it allows the ordering of the $n$ numerical values $x^{i}>0$ in the vector $x$ to become unambiguous by properly ordering the $n$ binary values $u^{i}_{\bar{k}}$ in the binary vector $u_{\bar{k}}$ for each value of $\bar{k}=1,...,T-1$.

The ambiguity in the ordering of the $n$ numerical values $x^{i}>0$ in $x$ given in (3) is implicitly included due to the fact

[TABLE]

and the fact that $x^{i}>0$. In this way, the values of $x_{i}$ can be ordered if we order the values of $u^{i}_{\bar{k}}\in[0,1]$ contained in the time ordered binary switch state vector $u_{\bar{k}}$ by the choice of a positive real linear function

[TABLE]

that satisfies the property

[TABLE]

The rationale behind the choice of the function $f(\cdot)$ in (15) is as follows. For a given value of $n$, there are $2^{n}-1$ possible binary combinations of the vector $x_{\bar{k}}$, excluding the value 0. Ordering the values of the time ordered binary switch state vector $u_{\bar{k}}$ according to (16), allows us to choose a fixed $x$ (independent of $\bar{k}$) with $x^{i}>0$ to satisfy (14).

An obvious choice for the desired function $f(\cdot)$ in (15) is to use the fact that $u^{i}_{\bar{k}}\in[0,1]$ and that the vectors $u_{\bar{k}}$ with length $n$ can be seen as a $n$ bit binary number representing a signed integer number $d_{\bar{k}}>0$. The conversion from an $n$ bit binary number $u_{\bar{k}}$ to a signed integer number $d_{\bar{k}}$ is given by

[TABLE]

which clearly satisfies the conditions of the function $f(\cdot)$ given in (15) and (16). With the choice of $f(u_{\bar{k}})=d_{\bar{k}}$ and ordering the integer numbers $d_{\bar{k}}$ according to (15), the static load optimization problem in (2) can be rewritten as

[TABLE]

with $f(u_{\bar{k}})$ given in (17). It should be noted that for a given value of $n$, the binary numbers and the ordering of $u_{\bar{k}}$ in (18) are completely known. As a result, only an optimization over $x$ is required reducing the optimization in (2) to an equivalent standard least squares (LS) optimization given in (18). With the known time ordered binary switch state vector $u_{\bar{k}}$, the time ordered error $E(t_{\bar{k}})$ for $\bar{k}=1,...,T$ in (2) can be written in a matrix notation $E=S-Ux$ where the matrices are given by

[TABLE]

and where the block rows of the matrix $U$ will be repeated entries given of the time ordered binary load switch vector \left[\begin{array}[]{cccc}u_{\bar{k}}^{1}&u_{\bar{k}}^{2}{}&\cdots&u_{\bar{k}}^{n}\end{array}\right] due to Kronecker product with the $1\times L$ unity vector $\mathbf{1}_{1\times L}$ where $L=T/(2^{n}-1)$. The repeated entries are needed to ensure that $U\in R^{T\times n}$, as there are only $2^{n}-1$ binary load switch combinations (excluding all loads off), while the available number of (solar) power data points $T>>2^{n}-1$. The typical value of $L$ for the repeating entries in the block rows of $U$ is $L=20$, using $20\cdot(2^{n}-1)$ points of ordered (solar) power data for computation of the optimal load size distribution. The standard LS minimization of (18) can be rewritten as

[TABLE]

where the solution can be computed by $\hat{x}=[U^{T}U]^{-1}[U^{T}S]$.

The standard LS solution in (18) will solve the Least Squares error of the static power mismatch, but does not ensure yet that the sum of the load distribution is bounded as in (4) to avoid oversizing of the loads in trying to match the anticipated maximum power production $S_{max}$ [43]. The additional linear equality constraint on the sum of the load distribution can easily be incorporated via an Equality Constrained Least Square (ECLS) problem

[TABLE]

that also includes a Lagrange multiplier $\lambda$ and the unit equality constraint vector $D=[1~{}1~{}\cdots~{}1]=\mathbf{1}_{1\times n}$ and a chosen value of $C$ in the range $0.5\leq C\leq 1$ to satisfy (4). The ECLS solution is now given by

[TABLE]

with $U$ and $S$ as given in (19). Since the optimal value $C$ in the range $0.5\leq C\leq 1$ to avoid oversizing of the loads by bounding the sum of the load distribution as in (4) is unknown, an additional line search along $C$ can be used to determine the optimal load oversizing constraint.

4.3 Inequality Constrained Least Squares Optimization

Although the ECLS approach presented above computes the optimal load size distribution $x=[x^{1}~{}x^{2}~{}\cdots~{}x^{n}]^{T}$, the optimal solution $\hat{x}$ to (20) still depends on the choice of the time ordered binary load switch vector $u_{\bar{k}}$ due to the bi-linear nature of the optimization problem in (2). Even when the (solar) power data $S(t_{k})$ has been ordered, it is still not clear when exactly the loads will be turned on/off. Referring to Figure 9 to illustrate this concept, it is not clear what the optimal load switch time $m_{1}$ for the first load $1$ will be and how many samples $m_{2}$ the first load should remain on before the second load is switched on at $m_{1}+m_{2}$ samples. Clearly, the optimal static values $\hat{x}^{1}$ and $\hat{x}^{2}$ of the loads depend on these switching times.

To explicitly incorporate the ordering constraint (3) and allow for variability in the switch time, both the optimization variable $x$ and the matrix $U$ in (19) can be modified., while still allowing for a convex optimization. First we define

[TABLE]

where the upper diagonal matrix $T$ ensures that $\bar{x}$ now reflect the incremental change in the load size distribution. In this way the ordering constraint (3) can be enforced explicitly by $\bar{x}>0$. Secondly, we define

[TABLE]

where now the block rows of $U$ are defined by

[TABLE]

in which the integer vector

[TABLE]

is the set of $2^{n}-1$ possibilities of incremental switching time values in the case of $n$ loads.

The incremental switching time values $m$ allow variable timing when loads are switched, similar as in Figure 9 for the case of $n=2$ loads. Moreover, given the vector of incremental switching time values $m$, the optimal incremental load size distribution can be computed with a Inequality Constrained Least Squares (ICLS) problem

[TABLE]

where the matrix $A$ and $b$ can be used to enforce inequality constraints. In particular, the choice

[TABLE]

enforces the ordering constraint (13) via $\bar{x}\geq 0$ and ensures load power demand is always under the (solar) power curve via $U(m)T\bar{x}=U(m)x\leq S$ to avoid oversizing of the loads directly.

The solution to the ICLS problem can be solved with standard convex optimization tools and will lead directly to optimal results for the static load distribution, given the integer vector

[TABLE]

of $2^{n}-1$ possibilities of incremental switching time values. An additional line search or non-linear optimization can be used on top of the ICLS problem to compute an optimal set of incremental switching time values to further improve the solar utilization of the static load distribution. Theoretically, such an additional search along the switching time values via an additional iterative or gradient based optimization should lead to the globally optimal solution to the bi-linear optimization problem of (2), as one can use the full number of $T$ data points on the (solar) power data, while using the smallest number of optimization variables. A drawback is that the iterative search for switching time values may get stuck in a local minimum. This was performed based on gradient search over $m$ using the Nonlinear programming solver (fmincon) in the MATLAB optimization toolbox

4.4 Mixed-integer Linear Programming (MILP) Method

Since the optimization problem (2) is a bi-linear optimization problem involving a mix of binary and real numbers, a limited range of optimization solvers can be applied and none of which may guarantee global optimality. To guarantee the optimality, we employ the Big-M relaxation method [44] to convert the optimization problem into a mixed-integer linear programming (MIPL) problem. The convexity of MILP therefore fulfills the zero duality gap. In addition, MILP has the capability of determining the exact switching schedule of load units, which the other approaches discussed in this paper lack. In fact, solving the optimization problem through MILP is equivalent to simultaneous solution to planning and scheduling problems.

There exist many mature MILP solvers which are capable of solving large-scale MILP problems with millions of variables within a reasonable time frame [45]. Proper selection of optimization variables and use of the disjunctive methods discussed in [22] make it possible to reformulate the original problem as an MILP problem.

Denoting the binary variable $u_{i}(t_{k})$ as the on/off status of the unit $i$ at time step $t_{k}$, the optimization problem (2) is presented as below,

[TABLE]

where $i$ and $t_{k}$ are the index for units and time steps respectively. $S(t_{k})$ denotes the solar power timeseries of length $T$. $x\in R^{n}$ is the vector of load sizes. $y_{i}(t_{k})$ also denotes the committed demand by load $i$ at time $t_{k}$.

The last constraint in (22) is a bi-linear equality constraint and makes (22) nonconvex. To resolve the nonconvexity issue, the big-M technique is used to replace the constraint with the two following sets of constraints,

[TABLE]

where $M$ is a big-enough positive number, e.g. $10^{6}$ as utilized in the numerical examples here. The constraint sets (23a) and (23b) are binding and relaxed respectively when $u_{i}(t_{k})=1$, guaranteeing $y_{i}(t_{k})={x_{i}}$. Likewise, $y_{i}(t_{k})$ exactly equal to zero $u_{i}(t_{k})=0$.

One great advantage of this method is that the input data can be one certain day (clear or cloudy), that is the data does not need to be sorted as compared to the other approaches. In addition, as no change is required in the order of PV power profile while solving the optimization problem via MILP, additional constraints and technologies such as minimum uptime and downtime of load units and employing energy storage systems could be considered in the optimization process. Since this paper is focused on optimal load sizing, these options will be discussed in detail in the future research works.

MILP was implemented using the CVX toolbox and Gurobi 6.50 to solve the integer problem in our optimization. To reduce computational expense, down-sampling the original data was needed. As the sequence of PV data does not affect the optimization results, down-sampling with the ratio of $1:n$ can be performed by arbitrarily selecting one data point from every $n$ data points from either the original or sorted data set further discussed in Section 5.3.

5 Results

In this paper, two levels of optimization are performed. The first level is toward unit sizing for a given number of units. After the size of the units is determined, they are scheduled through an optimization process to capture the maximum solar power. For example Figure 10, which will be discussed more in section 6, illustrates the scheduling results on three sample days for different number of units.

San Diego solar power data, described earlier in Section 3 is investigated in the following subsection using the different optimization techniques. This is followed by an additional case study and concluded with a discussion.

5.1 Equality Constrained Least Squares Optimization (ECLS) Results

The optimal unit sizes $x$ are obtained from equation (20). Subsequently, the operation is simulated for one year of solar power data. These results do not constrain the units with any minimum up or down time meaning the units are instantly turned on or off. After that the solar utilization ($SU$) was calculated by dividing energy consumed by the loads over the solar power data for one year.

A Sensitivity analysis for ECLS method was performed for $n=3$ showing the 75th percentile of the maximum solar utilization in 42 points and showing the range of the unit sizes. The maximum solar utilization for $n=3$ as shown in Table 1 is 0.8277 but for other $x$ combinations the solar utilization can be as low as 0.7542 (Figure 11). Unit size appears to be closely related to solar utilization. For example the size of $x_{1}$ varies between 0.5752 and minimum is 0.4385 with quartiles of 0.5418, 0.5068, 0.5418. This results in a 9% change in solar utilization with this unit size range. The smallest unit sizes are associated with the smallest solar utilization. Then there appears a bifurcation where the largest and mid-size units achieve medium solar utilization. The largest solar utilization are associated with medium-to-large unit sizes.

5.2 Inequality Constrained Least Squares (ICLS) Results

The Inequality Constrained Least Squares (ICLS) optimization results are shown in Table 2 and will be discussed in section 6.

5.3 MILP Optimization Results

To mitigate the adverse effect of down-sampling on the simulation results in this paper, down-sampled data points are selected uniformly from sorted data set, which represents the original data set more accurately. Figure 12 examines the impact of downsampling on computational speed and solar utilization for the case of 3 units $(n=3)$. The final solar utilization is the nearly the same for different down-sampling rates (with in 1% variation) while the computational cost varies remarkably. Table 3 shows the best solar utilization obtained in Figure 12.

For comparison Thuwal, Saudi Arabia case was performed using the ICLS approach and the results are shown in Figure 13. The resulting unit sizes are expected to be larger than San Diego due to the lack of cloudy days driving the need for smaller, more adaptive units.

6 Discussion

Looking at the solar utilization comparison between the three optimization approaches shows close results for each $n$ number of units as shown in Figure 14. Figure 14 shows summary of optimal unit size achieved by each method as well as solar utilization. ECLS was chosen as the reference case and assigned a solar utilization of 1 as shown in Figure 15. Theoretically the ICLS optimization approach should lead to the largest solar utilization, as it uses the full data set and the smallest number of optimization variables. ICLS indeed yielded the best results, but improvements compared to ECLS only ranged from 3.9% to 0.1% decreases as $n$ increases. However, MILP results were slightly less than ICLS. From Figure 14, it is also clear that as $n$ increases the different among the approaches decreases.

Figure 13 summarizes the differences between the unit sizes for both case studies (San Diego, USA and Thuwal, Saudi Arabia) as well as the optimal solar utilization obtained. Keep in mind that these results are normalized for both solar power and unit size. The main difference between the two sites is larger units are preferred in the Thuwal case with an average size increase of 25%. For $n=2$ the solar utilization for Thuwal was larger compared to San Diego, mainly due to a higher clear day count over the year. As $n$ increases both site solar utilization tend to get closer until they nearly match for $N=6$.

Each method proposed in this paper carries advantages and drawbacks when solving the unit sizing problem, as summarized in Table 4. From a global optimality point of view, the analytic approach is proven to achieve the global optimum, but it can not be generalized for all solar day patterns. The other approaches on the other hand were proven to be robust and suitable for real sizing problems. MILP approach guarantees global optimality, but its computation expense increases exponentially with the number of decision variables and constraints. ECLS required downsampling as well, and that explains the decreases in solar utilization difference between different methods as $n$ increases, where downsampling get more accurate. Finally, ICLS is not guaranteed to converge to the global optimal, rather could get stuck in a local minimum. Thus, the final solution highly depends on initial conditions. This could be avoided by slightly perturbing the initial or final solution and restarting the optimization and thereby exploring the solution space

The analytic approach can not be generalized for all solar day patterns; it requires a functional form for the solar data and function inverse has to be known which limits this approach to be applied to run yearly data. For that reason, the analytic method can not be used for real sizing problems. The MILP technique is the most flexible approach especially given the ability of adding more constraints to the problem, such as battery or minimum up and downtime. It can also solve a daily pattern which is suitable for schedule peruses with global optimality claims in case of convergence, still the number of variables play a big role in convergence and determining optimality. The main obstacles with the MILP are the high computational cost due to a large number of variables as well as the approximation or the relaxation applied to the problem. ECLS and ICLS are both scalable to solve multiple number of years with fast computation time. The main disadvantages of these approaches are that they do not consider the daily solar profile in addition to that they have limited capability to add constraints or improve the case as in adding batteries or minimum up and downtime.

Sample results were selected in Figure 10 for $n=2$ up to $n=6$ using the ICLS approach as shown in Table 2 and Figure 14. The algorithm solves for the optimal unit sizing over one year which contains various daily patterns based on the prevailing weather conditions. For the clear day (D1), clearly a large area of the peak of the day is wasted but lost energy is reduced as the number of units $n$ increases. On the other hand, on the most cloudy day D3, the loss was reduced. This is a reflection of the input data where clear midday periods that yield normalized power output close to 1 are less common than cloudy days and morning and evening periods. Consequently, the unit sizes are selected to track lower power outputs more closely and clear midday periods are curtailed. If the algorithm was intended to solve the optimization over just one of the days pattern or for a different site, the results would be different.

Moreover, the computational time for the different methods was performed in a 3.4 GHz Intel Core i7 processor with 32 GB of RAM. As discussed earlier the analytic approach can not solve the planning problem, but it is very fast routine since it only computes derivative of existing function and substitute variables. MILP is the slowest, for the full set of variables, whereas Figure 12 shows it takes 7000 seconds to solve for input with 1:50 samples ratio and only 170 second for 1:200 downsample ratio. As discussed before the ICLS optimization is expected to return the optimal sizing of units, it does not require removing the zero solar radiation data. The main advantage of ICLS is that it searches for the optimal unit size and switching time by avoiding the equidistant constraints giving by the ECLS method. The main disadvantage for ICLS, it is sensitive to the initial condition, especially for larger n ($n>4$) which means that global optimality is not guaranteed. The computation time for ICLS varies from 11 seconds for $n=2$ to around 150 seconds for $n=6$. ECLS is the fastest approach even though it has to search for the optimal C given in (20). ECLS forces the results of the switch time to be equidistant outdistance spaced from each other which is not the optimal decision. Also it requires removal of zeros in the input data; otherwise the results and the solar utilization will be negatively affected. Each run of ECLS costs 0.01 seconds so solution time is dependent on the resolution of the line-search (C). For our purposes, ECLS is the fastest approach. Assuming the resolution of the line-search is 100 steps, ECLS is 10 times faster than ICLS for $n=2$ and 140 faster for $n=6$.

7 Conclusions

Solar PV is a desirable energy source for many standalone applications. Variability of solar irradiance is one of the main challenges of PV utilization. State-of-the-art optimization techniques were developed and applied to optimize the solar utilization by sizing a given number of load units based on one year of data collected in San Diego and Saudi Arabia. The algorithm switches the units on and off to “load follow” the available solar power during the day to maximize the solar energy utilization. The primary output of the algorithms is the optimum sizing for a given number of units, but unit scheduling is a byproduct of the analysis. Three different optimization methods are proposed to solve for the optimal unit size: Equality Constrained Least Squares (ECLS), Inequality Constrained Least Squares (ICLS), and Mixed-Integer Linear Programming (MILP). The performance of the three methods was compared with two case studies. Results for the San Diego case indicate a solar utilization (i.e., percentage of energy captured by units over available energy) differed by less than 5% between the algorithms. It was shown the utilization increased from 73% for two units up to 98% for six units for San Diego case. As expected, the ICLS optimization yields the largest utilization. The results obtained will differ by location and may even vary year-to-year due to spatio-temporal patterns in the solar resources and cloud coverage as evident from the differing results between San Diego, USA and Thuwal, KSA. The methodology proposed in this paper allows computationally efficient solutions even when several years of solar resource data are available and yield the optimal sizing for the given data. For practical applications, the economics also need to be considered as smaller units typically cost more per kW and an optimization based on cost would therefore yield larger and prefer fewer units. Within our framework, it is possible to assign a cost function to the number of units and to the solar utilization to provide solutions for practical applications.