Unconstrained Estimation of Multitype Car Rental Demand

Abstract

The unconstrained demand forecast for car rentals has become a difficult problem for revenue management due to the need to cope with a variety of rental vehicles, the strong subjective desires and requests of customers, and the high probability of upgrading and downgrading circumstances. The unconstrained demand forecast mainly includes repairing of constrained historical demand and forecasting of future demand. In this work, a new methodology is developed based on multiple discrete choice models to obtain customer choice preference probabilities and improve a known spill model, including a repair process of the unconstrained demand. In addition, the linear Holt–Winters model and the nonlinear backpropagation neural network are combined to predict future demand and avoid excessive errors caused by a single method. In a case study, we take advantage of a stated preference and a revealed preference survey and use the variable precision rough set to obtain factors and weights that affect customer choices. In this case study and based on a numerical example, three forecasting methods are compared to determine the car rental demand of the next time cycle. The comparison with real demand verifies the feasibility and effectiveness of the hybrid forecasting model with a resulting average error of only 3.06%.

1. Introduction

The car rental industry plays a huge role globally. It can act as a lubricant for production and consumption to ease their mutual restraints. It can also expand the automotive consumer market and evaluate the popularity of new cars before they come into the market. At the same time, because the public transportation system is limited by operating time, departure frequency, accessible range, comfort, privacy, and other conditions, car rental, with its outstanding advantages such as high flexibility, ease of use, and privacy, has shown considerable growth over the years. The global car rental market was valued at $79.5 bn in 2018 and is expected to reach $125.4 bn by 2025, registering a CAGR of 5.1% from 2019 to 2025 [1]. In recent years, the demand for car rental services in the Asia-Pacific region has become the largest and fastest growing segment, especially in China and India, where there are high population density and rapidly growing demand.

Mainstream car rental companies started to use revenue management (RM) as a crucial tool for reducing operating costs and improving competitiveness in the 1990s. RM is an important interface between enterprises’ supply and market demand, and its objective is to maximize profits. The most widely known RM application occurs in the airline industry. It mainly includes four aspects, namely forecasting, overbooking, capacity control, and pricing [2]. Among them, forecasting is the foundation and provides vital input information for other work. A 12.5–25% reduction of forecast error could translate into 1–3% incremental revenue generated from the RM system [3,4]. Therefore, to achieve the desired level of accuracy, it is important and necessary to innovate the forecast method constantly [5]. In reality, influenced by booking and/or capacity limits, the customer’s present demand may not be satisfied, which will lead to it being lost (spilled) or recaptured by a more expensive (buyup) or cheaper (buydown) available product. The sample booking curves with three grades of products (the fare order is A > B > C) are shown in Figure 1 [6,7]. Therefore, in order to repair the constrained demand, it is necessary to carry out an unconstrained estimation of the customer’s actual demand and forecast the future demand on this basis [8,9].

Research into unconstrained estimation (sometimes also called detruncation, uncensoring) in the RM system started more than 40 years ago and mainly focused on the airline industry. Unconstrained estimation mainly uses observed “censored demand” to estimate and calculate the “uncensored demand”. The original methods were based on statistics, such as the expectation maximization (EM), projection detruncation (PD), booking curve (BC), and mean imputation. Since the beginning of this century, discrete customer choice models have grown in popularity to estimate unconstrained demand [8,10,11]. Forecasting usually yields the demands and no-show rate by booking class, the latter of which will not be discussed here. The common methods include time series, linear regression, exponential smoothing, double exponential smoothing, pickup, neural networks, principal component analysis, and adaptive models [6,12,13]. Since air transportation can meet the travel demand between two points, the prediction of OD demand is necessary [14,15,16]. Guo surveyed the history of research on unconstraining methods by reviewing over 130 references [17]. Weatherford provided two reviews of unconstrained estimation and forecast methods broadly used in different industries up to 2016 [18,19]. Recent discussions include forecast multipliers and hybrid forecasting [9], the effect of customers’ reference price on demand [20], Gaussian processes for unconstraining demand [21], and demand forecast accuracy [22]. Due to the important research and application value of revenue management theory, many advances have been made in the hotel, car or truck rental, cargo/freight, internet service and retailing, cruise line, rail, and other industries regarding demand estimation and forecasting [19].

Unlike the airline industry, the car rental industry is more concerned with the demand at each rental station, and OD forecasting does not need to be considered as much. Supply and demand are dynamically balanced, but there is still an overestimation or underestimation. Since demand determines the size of the leasing station inventory and fleet, demand estimation and forecasting techniques are worth studying [23]. The first reports of RM application to car rental were at Hertz and National Car Rental [24,25]. Jishan used a decomposition method to identify the actual latent demand from the total recorded turndowns [26]. Gordon simplified the problem by human–computer interaction, hoping to find the right level of detail at which to forecast and optimize [27]. Dhruval optimized the fleet with the help of a robust design methodology, considering fleet as a product and cars as an affecting variable [28]. Kourentzes addressed the frequently encountered situation of observing only a few sales events at the individual product level and proposed variants of small demand forecasting methods to be used for unconstraining [29]. In general, there is relatively little research on unconstrained estimation and prediction of car rental demand. By referring to research methods in other fields, relevant research in the car rental industry is worthy of a further discussion.

Due to the relative lack of articles that deal with car rental demand forecasting, this work develops an optimization methodology to estimate and forecast the demand for a car rental station; it is then applied to a case study in China. By using the case study, the work also enables comparison between the effects of different methods, especially for cases of small demand. The remainder of this work is organized as follows. Section 2 presents the problem statement and hypothesis, including an improved spill model for unconstrained estimation. In Section 3, we propose a hybrid forecasting method based on the Holt–Winters model and backpropagation neural network. Subsequently, Section 4 introduces the case study with numerical experiments and a discussion. Section 5 concludes the work.

2. Unconstrained Estimation

2.1. Hypothesis and Variable Definitions

Before studying the forecasting method, the conditional hypothesis and variable definitions should be explained. This work is mainly aimed at the car rental industry, and the method is also applicable to the hotel and aviation industries, etc.

Hypothesis 1.

Unconstrained demand $X$ of car type $i$ in presell lead time meets normal distribution, and the lead time is seven days.

Hypothesis 2.

Demands of various car types are random and independent, and the numbers of car rentals are counted by day; the deadline is midnight.

Hypothesis 3.

Inventory quantity is constant without car dispatching among stations, taking no account of batch demand, overdue return, or oversold behaviors.

Hypothesis 4.

The customer’s willingness to pay is from low to high.

The variables are defined as follows:

• $i$: index of car types, $i\in I$. Indexes of $I$ are listed in an increasing order of car rental price, corresponding to the gradually increased levels;
• $j$: index of rental demand recording times, $j=1,2,\cdots ,\mathrm{J}$;
• $t$: review points of inventory control, $t=0,1,\cdots ,\mathrm{T}$; presale system will open when $t=0$, and it will be closed when $t=\mathrm{T}$.
• ${\Delta }_t$: presale lead time interval between $t$ and $t-1$;
• $B(i,j,t)$: order quantity limit of car type $i$ in ${\Delta }_t$ in the $j\mathrm{th}$ documented circle, namely the quota ceiling for rent;
• $C(i,j,t)$: cumulative order quantity of car type $i$ at review point $t$ in the $j\mathrm{th}$ documented circle, particularly $C(i,j,0)=0$;
• $I(i,j,t)$: presale state of car type $i$ in ${\Delta }_t$ in the $j\mathrm{th}$ documented circle, “1” indicates that demand is not restricted, and “0” indicates that demand is limited.
• $IB(i,j,t)$: order quantity of car type $i$ in ${\Delta }_t$ in the $j\mathrm{th}$ documented circle, particularly $IB(i,j,0)=0$, $IB(i,j,t)=C(i,j,t)-C(i,j,t-1)$;
• $U(i,j,t)$: spillage of car type $i$ in ${\Delta }_t$ in the $j\mathrm{th}$ documented circle, particularly $U(i,j,0)=0$;
• $CP(i,j,t)$: complete spillage of car type $i$ in ${\Delta }_t$ in the $j\mathrm{th}$ documented circle, particularly $CP(i,j,0)=0$;
• $SU(i,j,t)$: cumulative spillage of car type 1 to $i$ at review point $t$ in the $j\mathrm{th}$ documented circle, $SU(i,j,0)=0$;
• $IU(i,j,t)$: unconstrained estimation of car type $i$ in ${\Delta }_t$ in the $j\mathrm{th}$ documented circle, particularly $IU(i,j,0)=0$;
• $TU(i,j,t)$: real value of car type $i$ in ${\Delta }_t$ in the $j\mathrm{th}$ documented circle.

2.2. Unconstrained Estimation Methods

For the revenue management system, the starting point of analysis is some set of assumptions regarding an underlying stochastic or deterministic demand process [30]. In fact, when the real customer demand exceeds the preset order quantity limit, the recorded demand data in the presale system would be censored. Using these “cutoff demand data” for demand forecasting is bound to reduce the accuracy of forecasts and directly impact the revenue. Therefore, deducing the real historical customer demand from the recorded demand data is the first problem. Some common unconstrained estimation methods are as follows.

2.2.1. Expectation Maximization

EM is one of the most common methods to restore incomplete data by using the iterative algorithm in statistics. The EM algorithm is mainly composed of two steps, E (expectation) and M (maximization)—that is, after initializing the parameters to be estimated through the observable incomplete data information, the conditional expectation of missing information in the incomplete data is calculated by parameter estimation, and the missing information is replaced accordingly. In step E, an expected log-likelihood function of complete data is obtained. In step M, the improved parameter values are obtained by using this function through the maximum likelihood estimation method. The whole process will be repeated until convergence [31].

Considering the unconstrained demand $X$ in $\Delta t$ with the expectation $\mu$ and variance ${\sigma }^2$, the initial values are then

$$\mu {(i,t)}^{(0)}=\frac{{\displaystyle \sum _{j=1}^{\mathrm{J}}IB(i,j,t)I(i,j,t)}}{{\displaystyle \sum _{j=1}^{\mathrm{J}}I(i,j,t)}}$$

(1)

$${\sigma }^2{(i,t)}^{(0)}=\frac{{\displaystyle \sum _{t=1}^{\mathrm{J}}{\left\{[IB(i,j,t)-\mu {(i,t)}^{(0)}]I(i,j,t)\right\}}^2}}{{\displaystyle \sum _{t=1}^{\mathrm{J}}I(i,j,t)-1}}$$

(2)

When all data is constrained, $\mu (i,t)$ and ${\sigma }^2(i,t)$ cannot be initialized properly, and the repairing process is abandoned. Let the unconstrained preliminary estimation value be equal to the increment of the observed value; then, the unconstrained estimation value after $k$ iterations will be

$$IU{(i,j_j,t)}^{(k)}=IB(i,j,t)+IU{(i,j,t)}^{(k-1)}$$

(3)

When only one set of data is constrained, ${\sigma }^2(i,t)$ initialization is abandoned, and the convergence is tested directly by $\mu (i,t)$. The $k\mathrm{th}$ iteration is divided into two steps.

• Step E: unconstrained estimation value is

  $$IU{(i,j,t)}^{(k)}=\mathrm{E}[X|X\ge IU{(i,j,t)}^{(k-1)}]=\frac{{\displaystyle {\int }_{IU{(i,j,t)}^{(k-1)}}^{\infty }xf(x)\mathrm{d}x}}{{\displaystyle {\int }_{IU{(i,j,t)}^{(k-1)}}^{\infty }f(x)\mathrm{d}x}}$$

  (4)

  where $f(x)$ is the normal probability density function, and $IU{(i,j,t)}^{(0)}=IB(i,j,t)$ when $k=0$.
• Step M: find the partial derivative of the log-likelihood function $Q(\mu ,{\sigma }^2)$ versus $\mu {(i,t)}^{(k-1)}$ and ${\sigma }^2{(i,t)}^{(k-1)}$, and let them equal 0 to calculate the re-estimated $\stackrel{⌢}{\mu }$ and ${\stackrel{⌢}{\sigma }}^2$.

  $$Q(\mu ,{\sigma }^2)=\ln L(\mu ,{\sigma }^2)=-\frac{\mathrm{J}}{2}\ln 2\pi -\mathrm{J}\ln \sigma -\frac{{\displaystyle \sum _{j=1}^{\mathrm{J}}{[IB(i,j,t)-\mu ]}^{2}I(i,j,t)}+{\displaystyle \sum _{j=1}^{\mathrm{J}}{[IB(i,j,t)-\mu ]}^2[1-I(i,j,t)]}}{2{\sigma }^2}$$

  (5)

  $$\stackrel{⌢}{\mu }=\mu {(i,t)}^{(k)}=\frac{1}{\mathrm{J}}\left\{{\displaystyle \sum _{j=1}^{\mathrm{J}}IB(i,j,t)I(i,j,t)}+{\displaystyle \sum _{j=1}^{\mathrm{J}}IU{(i,j,t)}^{(k-1)}[1-I(i,j,t)]}\right\}$$

  (6)

  $${\stackrel{⌢}{\sigma }}^2={\sigma }^2{(i,t)}^{(k)}=\frac{1}{\mathrm{J}-1}\{{\displaystyle \sum _{j=1}^{\mathrm{J}}{[IB(i,j,t)-\mu {(i,t)}^{(k-1)}]}^{2}I(i,j,t)}+{\displaystyle \sum _{j=1}^{\mathrm{J}}{[IU{(i,j,t)}^{(k-1)}-\mu {(i,t)}^{(k-1)}]}^2[1-I(i,j,t)]}\}$$

  (7)
  Convergence condition: the operation should be stopped when $|\mu {(i,t)}^{(k)}-\mu {(i,t)}^{(k-1)}|\le \epsilon$ or $|{\sigma }^2{(i,t)}^{(k)}-{\sigma }^2{(i,t)}^{(k-1)}|\le \epsilon$, where $\epsilon$ is any small given number. Otherwise, go back to step E and continue iterating.

2.2.2. Projection Detruncation

The initial values, convergence time, and iteration number of the EM algorithm in the iterative process are more susceptible to their own data in varying degrees. Based on the EM algorithm, Hopperstad proposed PD, which solves the limitation of EM in a large-scale unconstrained estimation [32,33]. The research proves that the PD algorithm is more flexible in practical applications.

PD and EM have the same principle, both of which include step E and M. Both algorithms use the iterative idea of statistics to estimate the parameters of constraint demand data, but the difference is that the parameter-related heuristic method is used in PD at step E. In PD, the unconstrained estimation value is gained without calculating the conditional expectation of the constraint data. In step M, the convergence conditions are the same, namely

$$\begin{array}{cc}\hfill IU{(i,j,t)}^{(k)}& =\mathrm{E}[X|X\ge IU{(i,j,t)}^{(k-1)}]\hfill \\ & ={\displaystyle {\int }_{IU{(i,j,t)}^{(k)}}^{\infty }f[x|\mu {(i,t)}^{(k-1)},{\sigma }^2{(i,t)}^{(k-1)}]\mathrm{d}x}\hfill \\ & =\tau {\displaystyle {\int }_{IB{(i,j,t)}^{(k-1)}}^{\infty }f(x|\mu {(i,t)}^{(k-1)},{\sigma }^2{(i,t)}^{(k-1)})\mathrm{d}x}\hfill \end{array}$$

(8)

When $IU{(i,j,t)}^{(k-1)}>\mu {(i,t)}^{(k-1)}$,

$$IU{(i,j,t)}^{(k)}=\mu {(i,t)}^{(k-1)}+\sqrt{-1.6{\sigma }^2{(i,t)}^{(k-1)}\ln \left\{1-{\left[1-2\tau +2\tau \Phi \left(\frac{IU{(i,j,t)}^{(k-1)}-\mu {(i,t)}^{(k-1)}}{\sigma {(i,t)}^{(k-1)}}\right)\right]}^2\right\}}$$

(9)

When $IU{(i,j,t)}^{(k-1)}\le \mu {(i,t)}^{(k-1)}$,

$$IU{(i,j,t)}^{(k)}=\mu {(i,t)}^{(k-1)}-\sqrt{-1.6{\sigma }^2{(i,t)}^{(k-1)}\ln \left\{1-{\left[1-2\tau +2\tau \Phi \left(\frac{IU{(i,j,t)}^{(k-1)}-\mu {(i,t)}^{(k-1)}}{\sigma {(i,t)}^{(k-1)}}\right)\right]}^2\right\}}$$

(10)

where $\Phi (x)$ is the distribution function of normal distribution, and parameter $\tau$ is any given constant ($0\le \tau \le 1$) and stays the same throughout the whole process. $\tau$ represents the probability that the observed constrained demand has a true value greater than $\mathrm{E}[X|X\ge IU{(i,j,t)}^{(k-1)}]$. Instead of calculating a conditional expectation, in the probability distribution graph (Figure 2), PD balances two things. The first is the area of the probability distribution between the order quantity limit (that is, the original constrained value) and the new estimate or projected value (area M), and the other is the area between the new estimate and infinity (area N). The horizontal axis value corresponding to region N represents the portion of the real demand that is underestimated by the unconstrained estimation [34]. When $\tau =0.5$, the PD algorithm can truly balance the values of M and N, yielding nearly similar results to the EM algorithm.

2.2.3. Multitype Spill Model

The multitype spill model (MSM) takes into account the vertical recapture of different price grade requirements, where the buyer’s behaviors of buyup and buydown indwell in different price levels. Therefore, the unconstrained estimation process based on MSM is closer to reality, which can effectively avoid the overestimation of real demand due to the repeated records of the revenue management system’s actual observable demand between different car type prices. The spill model is mainly divided into two parts, the first part being the calculation of spillage, and the second part being the calculation of unconstrained demand [35]. The specific calculation process is as follows:

• Check if the rental demand is constrained in ${\Delta }_t$. When $I(i,j,t)=1$, the rental demand is not constrained and no repair is performed; otherwise, the repairing process is performed, go to next step.

  $$I(i,j,t)=\{\begin{array}{c}1\hspace{0.17em},\hspace{0.17em}C(i,j,t)-B(i,j,t)<0\\ 0\hspace{0.17em},\hspace{0.17em}C(i,j,t)-B(i,j,t)\ge 0\end{array}$$

  (11)
• Parameters initialization. The distribution parameters of the demand in ${\Delta }_t$ are initialized by the observable “unconstrained demand”, which is the same as the EM algorithm.
• Calculation of spillage and cumulative spillage. According to the assumptions, the customer’s willingness to pay is characterized from low to high car type, so the spillage should be calculated from the lowest car type.
  • For a single car type, cumulative spillage can be determined by the following formula:

    $$SU(1,j,t)=\{\begin{array}{cc}\hfill 0\hfill & ,\hspace{0.17em}I(1,j,t)=1\hfill \\ \hfill {\displaystyle {\int }_c^{\infty }{f}_{i,t}(x)(x-c)\mathrm{d}x}\hfill & ,\hspace{0.17em}I(1,j,t)=0\hfill \end{array}$$

    (12)

    where $c=IB(i,j,t)$, $f_{i,t}(x)$ is the probability density function of unconstrained demand for car type $i$ in ${\Delta }_t$. When $I(1,j,t)=0$, rental demand was constrained,

    $$SU(1,j,t)={\displaystyle {\int }_c^{\infty }{f}_{i,t}(x)(x-c)\mathrm{d}x}=\sigma \varphi (b)+(\mu -c)[1-\hspace{0.17em}\Phi (b)]$$

    (13)

    where $b=\frac{c-\mu }{\sigma }$, $\varphi (x)$ is the distribution function of standard normal distribution, and $\Phi (x)$ is the distribution function of standard normal distribution [36].
  • For a multitype situation, when $i=1$, considering cumulative spillage $SU(I,j,t-1)$ at review point $t-1$ from car type 1 to $I$.

    $$SU(I,j,t)=\{\begin{array}{cc}\hfill SU(I,j,t-1)\hfill & , I(i,j,t)=1\hfill \\ {\displaystyle {\int }_c^{\infty }{f}_{1,t}(x)[x-c+SU(I,j,t-1)]\mathrm{d}x}& , I(i,j,t)=0\hfill \end{array}$$

    (14)
    When $I(i,T_j,1)=0$,

    $$SU(I,j,t)={\displaystyle {\int }_c^{\infty }{f}_{1,t}(x)[x-c+SU(I,j,t-1)]\mathrm{d}x}=\sigma \varphi (b)+[SU(I,j,t-1)-\sigma b][1-\Phi (b)]$$

    (15)
    When $i>1$, calculate $SU(i,j,t)$ by considering the cumulative spillage $SU(i-1,j,t)$ at review point $t$ from car type 1 to $i-1$.

    $$SU(i,j,t)=\{\begin{array}{cc}\hfill SU(i-1,j,t)\hfill & , I(i,j,t)=1\hfill \\ \hfill {\displaystyle {\int }_c^{\infty }{f}_{i,t}(x)[x-c+SU(i-1,j,t)]\mathrm{d}x}\hfill & , I(i,j,t)=0\hfill \end{array}$$

    (16)
    When $I(i,T_j,1)=0$,

    $$SU(i,j,t)={\displaystyle {\int }_c^{\infty }{f}_{i,t}(x)[x-c+SU(i-1,j,t)]\mathrm{d}x}=\sigma \varphi (b)+[SU(i-1,j,t)-\sigma b][1-\Phi (b)]$$

    (17)
• Calculate spillage.

  $$U(i,j,t)=\{\begin{array}{ll}SU(1,j,t),& i=1\\ SU(i,j,t)-SU(i-1,j,t),& i>1\end{array}$$

  (18)
• Get the repaired result, namely unconstrained estimation.

  $$IU(i,j,t)=\{\begin{array}{ll}IB(i,j,t),& I(i,j,t)=1\\ IB(i,j,t)+U(i,j,t),& I(i,j,t)=0\end{array}$$

  (19)

2.3. Improved Multitype Spill Model

Generally, traditional revenue management assumes that customer selection behavior is short sighted and ignores customer subjectivity. However, competition in the car rental market has intensified, and the “buyer market” has gradually expanded, so the subjective initiative of customers has been more fully reflected. The customer can choose target commodity according to his or her own preferences and subjective utility. These factors can be divided into three categories according to the differences in sources, namely accidental factors, objective factors, and subjective factors.

Accidental factors are not expected and have some randomness, such as concerts, sports events, and weather changes. Objective factors are external factors, which are generally not changeable in a short period of time. Objective factors are not directly related to the customers themselves, including policy and regulations, regional economic level, and rental vehicle attributes. Subjective factors are the customer’s own determinants, which are expressed in terms of customer preferences, trip purposes, age, and literacy. Subjective factors can be grasped by the analysis of the customer’s choice intention, usually in the form of a questionnaire.

2.3.1. SP/RP Survey

The basic data of customer selection preference were obtained by stated preference (SP) and revealed preference (RP) surveys. The SP survey preresearches the subjective preferences of customers when renting a car before the rental occurs. The RP survey analyzes the actual selection behavior of customers in the context of what has happened. SP data were more flexible, while RP data were more reliable.

The SP survey can obtain multiple data from a single interviewee, with the advantages of small sample and low survey cost. The basic principle of the SP survey is to predetermine the various attribute factors and their influence levels and then invite interviewees to score the plan in the set situation. The respondents’ choices may be different from their actual actions, but they can still find the main factors that can be calculated and have a large impact on the overall mean. Then, the secondary factors that are difficult to measure will be removed. The data obtained from the RP survey are the actual occurrence data, since the data records behaviors that have been taken by the respondent or the actual selection behavior observed. The basic principle of the RP survey is to describe the occurrence or existence of the situation with different attribute factors and influence levels and to invite the respondents to score the content.

Before the SP survey plan is designed, condition attributes, decision attributes, and corresponding level values should be determined. The survey content includes the importance of the brand effect, the satisfaction of car rental type, the sensitivity of the rental price, returning location, satisfaction of the vehicle status, and the satisfaction with a special offer, service attitude, procedure convenience, and accident handling. The RP survey is aimed at customers who have made choices, and the influencing factors are rated as very satisfied, satisfied, general, not satisfied, and very dissatisfied. The survey includes the type of rental vehicle and satisfaction with the rental price, return location restrictions, vehicle status, service attitude, procedure convenience, and accident handling when the preferred rental vehicle is rejected. See Appendix A and Appendix B for the SP/RP questionnaire of the case study.

Customers are affected by many attributes when selecting multiple models. In order to determine the degree of influence of each attribute, AHP is usually used. AHP (analytic hierarchy process) is an easy way to quantitatively and qualitatively deal with some fuzzy and complex subjective problems. The most critical step is to obtain a judgment matrix by pairwise comparison to calculate the weight of each factor. The specific calculation process is not given. However, due to the uncertainty of expert experience, the subjective constructed judgment matrix may not meet the conditions of complete consistency. In order to reduce the inconsistency of subjective judgment, the VPRS (variable precision rough set) can be used to simplify the factors affecting the customer’s choice behavior, and similar factors are described by several factors to construct the judgment matrix.

VPRS expands the Pawlak rough set by introducing a confidence level $\eta \left(0.5<\eta \le 1\right)$, so that admissible classification error is allowed to a certain extent. It can improve the flexibility of decision rules and, at the same time, reflect the correlation in data analysis, which is beneficial to the discovery of relevant data from unrelated data, that is, the implicit patterns in the data can be more clearly expressed [37].

Let $S=(U,A,V,f)$ be an information system, where $U=\{u_1,u_2,\cdots ,u_{|U|}\}$ is a finite nonempty set called the domain object space. $A=\{a_1,a_2,\cdots ,a_{|A|}\}$ is a finite nonempty subset of attributes. If the attribute in $A$ can be further divided into two disjoint nonempty subsets, namely the conditional attribute set $C\ne Ø$ and decision attribute set $D\ne Ø$, $A=C\cup D$, $X\subseteq U$, $C\cap D=Ø$, then $S$ is called the decision table. $V$ is the union of the attribute domains, i.e., $V=\cup V_a$ for $a\in A$, where $V_a$ denotes the domain of attribute $a$; $f$ is an information function which associates a unique value of each attribute with every object belonging to $U$.

For any subset $P\ne Ø$ and $P\subseteq C$, $R(P)=\{(x,y)\in U\times U:f(x,p)=f(y,p),\forall p\in P\}$ is an equivalence relation on $U$, then the corresponding equivalence class is represented as ${[x]}_p$. By replacing the inclusion relation with a majority inclusion relation in the original definition of lower approximation of a set, we obtain the following generalized notion $\eta$-lower approximation and $\eta$-upper approximation:

$${\underset{\_}{apr}}_p^{\eta }(X)=\cup \left\{{[x]}_p:\frac{\left|{[x]}_p\cap X\right|}{\left|{[x]}_p\right|}\ge \eta \right\}$$

(20)

$${\overline{apr}}_p^{\eta }(X)=\cup \left\{{[x]}_p:\frac{\left|{[x]}_p\cap X\right|}{\left|{[x]}_p\right|}\ge 1-\eta \right\}$$

(21)

Then, the classification accuracy can be defined as

$${\gamma }^{\eta }(P,D)=\frac{\left|\cup \left\{x\in U:\frac{\left|{[x]}_p\cap X\right|}{\left|{[x]}_p\right|}\ge \eta \right\}\right|}{\left|U\right|}$$

(22)

where ${\gamma }^{\eta }(P,D)$ measures the percentage of knowledge that can be correctly classified in the existing knowledge given a certain $\eta$ in the domain.

2.3.2. Intention Analysis

Nine factors affecting the customer’s car rental behavior are regarded as conditional attributes by the Delphi method (

Table 1. Factors affecting customer car rental choice behavior.

Factors	Subjective Factors	Objective Factors	Random Factors
Company brand $c_1$		√
Car rental type $c_2$	√
Rental price $c_3$	√
Special offer $c_4$			√
Service attitude $c_5$	√
Returning location $c_6$		√
Vehicle status $c_7$		√
Procedure convenience $c_8$		√
Accident handling $c_9$			√

), and the level values of the conditional attribute are defined by three grades: 1 = satisfaction, 2 = something in between, and 3 = dissatisfaction. The level values of the decision attribute can also be defined by three grades: R = rent, B = book but not necessarily rent, and G = give up.

Two rental sites (International Convention & Exhibition Center and Jiangbei Airport) of China Auto Rental Inc. (Beijing, China) were selected in Chongqing, China. Interviewees were randomly selected near each site, and the investigators cooperated to fill out the questionnaire. A total of 50 questionnaires were distributed and 40 questionnaires were returned, of which 30 were valid. By sorting out and analyzing the survey data, customer selection decision information can be obtained, as shown in

Table 2. Customer selection decision information.

Object	Condition Attribute									Decision Attribute
	${\mathit{c}}_1$	${\mathit{c}}_2$	${\mathit{c}}_3$	${\mathit{c}}_4$	${\mathit{c}}_5$	${\mathit{c}}_6$	${\mathit{c}}_7$	${\mathit{c}}_8$	${\mathit{c}}_9$
$u_1$	3	2	1	2	1	2	2	1	2	B
$u_2$	3	1	1	1	1	1	1	2	2	G
$u_3$	2	1	2	2	2	2	2	1	1	R
$u_4$	1	1	1	1	1	2	1	1	1	B
$u_5$	1	2	2	3	2	2	3	1	3	G
$u_6$	2	3	3	3	3	3	3	3	3	G
$u_7$	1	3	2	3	2	1	1	1	1	R
$u_8$	3	1	1	1	1	1	1	2	2	R
$u_9$	3	2	1	2	1	2	2	1	2	B
$u_{10}$	2	3	3	3	3	3	3	3	3	G
$u_{11}$	2	2	1	1	1	2	1	1	2	R
$u_{12}$	3	1	1	1	1	1	1	2	2	R
$u_{13}$	3	2	3	3	3	2	3	1	3	G
$u_{14}$	1	3	2	3	2	1	1	1	1	R
$u_{15}$	3	1	1	1	1	1	1	2	2	B
$u_{16}$	3	2	1	2	1	2	2	1	2	B
$u_{17}$	1	1	1	1	1	2	1	1	1	B
$u_{18}$	3	1	1	1	1	1	1	2	2	G
$u_{19}$	3	2	1	2	1	2	2	1	2	B
$u_{20}$	2	1	2	2	2	2	2	1	1	R
$u_{21}$	1	3	2	3	2	1	1	1	1	R
$u_{22}$	1	1	1	1	1	2	1	1	1	B
$u_{23}$	2	1	2	2	2	2	2	1	1	R
$u_{24}$	2	3	3	3	3	3	3	3	3	G
$u_{25}$	1	1	1	1	1	2	1	1	1	B
$u_{26}$	2	3	3	3	3	3	3	3	3	G
$u_{27}$	1	3	2	3	2	1	1	1	1	R
$u_{28}$	2	1	2	2	2	2	2	1	1	R
$u_{29}$	3	1	1	1	1	1	1	2	2	B
$u_{30}$	3	2	1	2	1	2	2	1	2	B

.

The following equivalence classes are available:

• $U/C=(X_1,X_2,\cdots ,X_9)$, where $X_1=\{u_1,u_9,u_{16},u_{19},u_{30}\}$, $X_2=\{u_2,u_8,u_{12},u_{15},u_{18},u_{29}\}$, $X_3=\{u_3,u_{20},u_{23},u_{28}\}$, $X_4=\{u_4,u_{17},u_{22},u_{25}\}$, $X_5=\left\{u_5\right\}$, $X_6=\{u_6,u_{10},u_{24},u_{26}\}$, $X_7=\{u_7,u_{14},u_{21},u_{27}\}$, $X_8=\left\{u_{11}\right\}$, $X_9=\left\{u_{13}\right\}$.
• $U/D=(D_{\mathrm{R}},D_{\mathrm{B}},D_{\mathrm{G}})$, where $D_{\mathrm{R}}=\{u_3,u_7,u_8,u_{11},u_{12},u_{14},u_{20},u_{21},u_{23},u_{27},u_{28}\}$, $D_{\mathrm{B}}=$ $\{u_2,u_5,u_6,u_{10},u_{13},u_{18},u_{24},u_{26}\}$, $D_{\mathrm{G}}=\{u_1,u_4,u_9,u_{15},u_{16},u_{17},u_{19},u_{22},u_{25},u_{29},u_{30}\}$.
• ${\underset{\_}{apr}}_p^{\eta }(D_{\mathrm{R}})=\{X_3,X_7,X_8\}$, ${\underset{\_}{apr}}_p^{\eta }(D_{\mathrm{B}})=\{X_1,X_4\}$, ${\underset{\_}{apr}}_p^{\eta }(D_{\mathrm{G}})=\{X_5,X_6,X_9\}$.

Classification accuracy of VPRS is an evaluation of the ability to perform object classification with a confidence level. Let $\eta =0.8$. Similarly, the classification quality of each attribute division can be determined according to Equations (20)–(22).

${\gamma }^{\eta }(P,D)=4/5$, ${\gamma }_{c_1}^{\eta }(P,D)={\gamma }_{c_2}^{\eta }(P,D)={\gamma }_{c_4}^{\eta }(P,D)=0$, ${\gamma }_{c_3}^{\eta }(P,D)={\gamma }_{c_5}^{\eta }(P,D)=7/15$, ${\gamma }_{c_6}^{\eta }(P,D)={\gamma }_{c_8}^{\eta }(P,D)=2/15$, ${\gamma }_{c_7}^{\eta }(P,D)={\gamma }_{c_9}^{\eta }(P,D)=1/5$.

The case where there is zero in the above classification quality indicates that some interleaved repetitions or substitutions exist among the nine main attributes. Therefore, there are six main factors affecting customer choice, namely rental price $c_3$, service attitude $c_5$, returning location $c_6$, vehicle status $c_7$, procedure convenience $c_8$, and accident handling $c_9$. By comparing the six elements in pairs, the importance of each factor in the conditional attribute on the customer decision can be obtained. The judgment matrix $\Lambda$ is

$$\Lambda =\left[\begin{array}{cccccc}1& 1& 7/2& 7/3& 7/2& 7/3\\ 1& 1& 7/2& 7/3& 7/2& 7/3\\ 2/7& 2/7& 1& 2/3& 1& 2/3\\ 3/7& 3/7& 3/2& 1& 3/2& 1\\ 2/7& 2/7& 1& 2/3& 1& 2/3\\ 3/7& 3/7& 3/2& 1& 3/2& 1\end{array}\right]$$

Obviously, the judgment matrix $\Lambda$ is a positive and negative matrix, which fully satisfies the consistency condition, and its maximum eigenvalue is ${\pi }_{\max }=4$. The corresponding feature vector is $\omega ={(0.629,0.629,0.180,0.269,0.180,0.269)}^{\mathrm{T}}$. After normalizing, the weight vector of $\Lambda$ is $\overline{\omega }={(0.292,0.292,0.083,0.125,0.083,0.125)}^{\mathrm{T}}$.

2.3.3. Customer Choice Probability

The multitype spill model is based on traditional revenue management. The assumption is that people in reality are irrational. When the product price is lower than the willingness to pay, the customer will buy it immediately. Obviously, this assumption ignores the initiative of the decision-making subject. In reality, individuals have various cognitive deviations. If we do not proceed from the customer’s choice behavior, we cannot find the customer’s car rental rules and the changing trend of demand, which is very unfavorable to the car rental company’s ability to set prices and inventory in the later stage.

Therefore, in order to describe the customer’s choice behavior, the Logit model is used to solve the above problems. The utility maximization theory is assumed, and the utility size is used to measure the probability that a customer selects a certain type of vehicle. The various attributes of different types constitute a comprehensive utility value. It is generally believed that the greater the utility of a customer to select a particular type of vehicle, the greater the probability that such a vehicle will be selected. Customers choose a variety of types and make a choice after evaluating the utility based on various attributes. The customer’s choice behavior can be described by a utility function $U_i=v_i+{\xi }_i$, where $v_i$ and ${\xi }_i$ are the average utility and the random utility error, respectively, corresponding to the customer’s selecting.

Assuming that the customer is rational, then each customer will choose the product that will maximize utility. The probability that the customer chooses car type $i$ is

$${\mathrm{P}}_i=\mathrm{P}(v_i\ge v_h,\forall h\in I,h\ne i)=\mathrm{P}(v_i+{\xi }_i\ge v_h+{\xi }_h,\forall h\in I,h\ne i)=\mathrm{P}(v_i+{\xi }_i\ge \underset{\forall h\in I,h\ne i}{\max }(v_h+{\xi }_h))$$

(23)

Assuming that ${\xi }_i$ are mutually independent and obey Gumbel distribution, then the probability variable deviations of the two independent Gumbel distributions are subject to the Gumbel distribution, which yields the general form of the multinomial logit model:

$$P_i=\frac{e^{v_i}}{{\displaystyle \sum _{i\in I}{e}^{v_i}}}$$

(24)

Since the RP survey is based on actual conditions, the survey respondents are customers who return cars to the rental sites. In the design of the questionnaire, the actual upgrade or downgrade has been fully reflected. Therefore, the corresponding utility when the customer’s rental car type is upgraded or downgraded is

$$U_{i,l}=0.292c_3+0.292c_5+0.083c_6+0.125c_7+0.083c_8+0.125c_9$$

(25)

If the customer did not lease type $i$ and then rented the type $h$, the probability of upgrading or downgrading is

$${\mathrm{P}}_{i,h}=\frac{e^{v_{i,h}}}{{\displaystyle \sum _{\begin{array}{l}l=0\\ l\ne h\end{array}}^I{e}^{v_{i,l}}}}$$

(26)

where ${\mathrm{P}}_{i,h}$ is the preference probability of buyup or buydown with car type $h$ when the customer’s choice of car type $i$ is rejected, and ${\mathrm{P}}_{i,0}$ is the leaving probability when car type $i$ is rejected.

The improved multitype spill model (IMSM) needs to calculate the complete spillage of all constrained data first and then calculate the demand transferring from other car types which are constrained according to the customer preference probability. Excluding this part is the spillage; then, the sum of cumulative order quantity and spillage will be the unconstrained estimation.

$$CP(i,j,t)=\{\begin{array}{ll}0,& I(i,j,t)=1\\ {\displaystyle {\int }_c^{\infty }{f}_{i,t}(x)(x-c)dx,}& {\displaystyle I(i,j,t)=0}\end{array}$$

(27)

Excluding the transferring is the spillage:

$$U(i,j,t)=CP(i,j,t)-[CP(1,j,t)P_{1,i}+CP(2,j,t)P_{2,i}+\cdots +CP(I,j,t)P_{I,i}]$$

(28)

Then, the unconstrained estimation is

$$IU(i,j,t)=IB(i,j,t)+U(i,j,t)$$

(29)

3. Hybrid Forecasting

3.1. Selection of Forecast Method

For demand forecasting, there are three methods: quantitative analysis, qualitative analysis, and decision analysis. Quantitative analysis relies on a large amount of historical data and can be divided into time series forecasting, which uses time to organize data, and the causal analysis method, which uses relationships to organize data. The time series method is the most commonly used, in which the Holt–Winters model and moving average method are representative. Time series forecasting can predict future demand, but it cannot explain the reason. Causal analysis commonly uses regression analysis and simulation methods. The causal analysis method uses the causal relationship between data to look for changes and is generally used for macro prediction.

Qualitative analysis mainly relies on expert knowledge and experience to evaluate and does not involve quantitative analysis. Delphi and judgments fall into this category. Although this kind of method operates simply, the subjectivity is too strong and the effect is also not favorable. The decision analysis method combines quantitative and qualitative methods. At present, the market survey and randomness method are more commonly used, but the research is not very mature. The quantitative analysis method is commonly used to compare the above three methods.

The Holt–Winters model is a kind of time series forecasting model that avoids the deficiencies of the moving average method. It uses a cubic smoothing equation to make different data have different weights, and the predicted value is the weighted sum of the previous data sequence. Larry Weatherford et al. deemed the availability of the basic neural network more useful than traditional forecasting methods (moving averages, exponential smoothing, linear regression, etc.) by comparing the mean absolute percentage error [38].

The backpropagation (BP) neural network is a widely used nonlinear forecast method that can simulate the neural structure of the human brain and solve more complicated problems. The BP neural network can use its own nonlinear characteristics to simulate the development trend of the data, without requiring an assumption function. When the prediction accuracy is reached, the future demand can be predicted according to the learning situation. However, the BP neural network relies on the initial conditions and is prone to fall into the local optimal solution. Therefore, the single prediction method has inevitable defects, and the error accuracy may not reach the conditions for actual use. This paper intends to use the hybrid forecasting to combine weights of multiple forecast methods and improve the forecasting accuracy.

3.2. Holt–Winters Model

The Holt–Winters forecasting model is observed to outperform other techniques for the time series, having changing seasonality, mean, and growth rate [39]. It is an adaptive model that automatically recognizes changes in data patterns. For example, if the deviation is caused by internal interference, it can be considered that the new observation has the same influence as the original data, and thus gives the same weight to the data of different periods. If the deviation is caused by external interference, then new observations and the original data have different influences on the prediction results, and the new observations have a higher impact on the prediction events. In order to show that the value of the data in different periods has a different influence on the forecasting results, different weights can be given to different periods.

The demand for car rental is a nonstationary time series with seasonal and cyclical trends. The Holt–Winters model can perform very accurate forecasting of this regular time series data, especially for trends and seasonal changes. It can decompose linear time series, seasonal variations, and random variation time series and properly filter the impact of random fluctuations. The Holt–Winters model is an improvement and development of the moving average method. It does not need to store much historical data but also considers the importance of each period of data and uses all historical data.

First, smoothing equations can be obtained by iteration:

$$\begin{array}{l}{T}_{i,t}=\alpha \frac{F_{i,t}}{U_{i,t-L}}+(1-\alpha )(T_{i,t-1}+S_{i,t-1})\\ S_{i,t}=\beta (T_{i,t}-T_{i,t-1})+(1-\beta )S_{i,t-1}\\ U_{i,t}=\gamma \frac{F_{i,t}}{T_{i,t}}+(1-\gamma )U_{i,t-L}\end{array}$$

(30)

where $T_{i,t}$ is the smoothed value, $S_{i,t}$ is the trend value, and $F_{i,t}$ is the seasonal index of car type $i$ at $t$ after iterations. $L$ is the length of the seasonal period, and $\alpha$,$\beta$,$\gamma$ are smoothing coefficients.

Then, the predicted value at the future time point $k$ is

$${\widehat{F}}_{i,t+k}=(T_{i,t-1}+k{S}_{i,t-1})U_{i,t-L+k}$$

(31)

where $U_{i,t-L+k}$ is the seasonal index of car type $i$, which is one cycle ahead of $t$, $k\in \{1,2,\cdots ,kL\}$.

The way to determine the three smoothing coefficients $\alpha$, $\beta$, and $\gamma$ is to minimize the error between the forecasted and actual values. In order to obtain more accurate and objective parameters, the traditional method is the residual square sum minimum method. Smoothing coefficients are all located in the interval (0, 1) and increase with 0.1 step length. The squares of the prediction residuals are calculated separately and summed until the smoothing coefficients corresponding to the sum of squares of the smallest residuals are found.

3.3. BP Neural Network

The neural network generally includes an input layer, an output layer, and a hidden layer. The input layer is located in the first, and no neurons are connected at the front end. The output layer and the hidden layer are all connected with neurons, and the weights have a one-to-one correspondence. The impact factors of output data include input, weight, threshold, and excitation function. Because of the large number of neurons, a large amount of information is stored, giving the neural network powerful data processing capabilities. The neural network has the advantages of a strong parallel processing ability, strong nonlinear processing ability, strong self-adaptation and learning ability, and strong associative memory and fault tolerance.

The BP neural network has a simple structure, strong plasticity, clear learning steps, and mathematical meaning. It has been proved that the BP neural network can simulate any complex nonlinear mapping by selecting three layers. The BP neural network is a feedforward type network that utilizes an algorithm of error back propagation. There are only feedforward associations between neurons, and no feedback, intralayer, or interlayer correlation. Linear sigmoid-type functions are generally used as excitation functions. Since the excitation function is measurable everywhere, for a BP network, the divided area is no longer a linear partition but an area composed of a nonlinear hyperplane, which is a relatively smooth surface. Compared with linear partitioning, this classification is more accurate and has greater fault tolerance. The learning method of the BP neural network is to strictly adopt the gradient descent method, so that the analytical formula of weight correction is also very clear.

3.3.1. Network Structure

A full connection is achieved between the upper and lower layers of the BP neural network, and there is no connection among each neuron. Studies have shown that when the output layer and input layer use a linear activation function and the hidden layer uses the sigmoid activation function, a BP neural network with a hidden layer can map all continuous functions. Therefore, when constructing a BP neural network model, only one hidden layer is generally used, as shown in Figure 3.

3.3.2. Learning Process

The learning process of the BP neural network consists of two parts: forward and reverse propagation. When the propagation is positive, the sample is delivered to the input layer and then processed by the hidden layer and passed to the output layer. If the value obtained by the output layer does not satisfy the expectation, then it enters the backpropagation link, which uses the error back to the input layer and continuously corrects the weight between layers during the transfer process. After repeated propagation, when the error is small enough to be acceptable, the learning process then stops (shown in Figure 4).

The specific learning steps are as follows:

Step 1

Program initialization. Select sigmoid as the activation function, then determine the minimum error, learning rate, and momentum coefficient.

Step 2

Calculate output. Input the initial weight and calculate the output values of the hidden layer and the output layer processing unit.

Step 3

Calculate the error value. When error is less than the given minimum error, go to Step 5; otherwise, go to Step 4.

Step 4

Backpropagation to adjust the weight between hidden layers; then reuse Step 2.

Step 5

Acquire the optimal output value and the weight of each layer and end the algorithm.

3.3.3. Sample Selection

After the model is established, the sample needs to be trained. In general, the training sample needs to meet four characteristics. (1) There is a close functional relationship between the input and object variables, and the object variable will change obviously due to the change of input variables. (2) Input variables are independent of each other. It is impossible to accurately calculate other components by using the components of input variables. (3) The data to be predicted has a certain commonality with the sample data. (4) Sample size should have a certain scale, so that the combination of all samples can reflect the mapping relationship between output variables and target variables.

The BP neural network applies the data-driven idea, that is, using a nonlinear characteristic to approximate a time series, and then using the clear logical relationship and historical data to express future values. Suppose there is a time series $\left\{X_i\right\}$, where the historical data is $X_p,X_{p+1},\cdots ,X_{p+q}$; to forecast the future value at time $p+q+r\hspace{0.17em}(r>0)$, some kind of nonlinear function relationship between $X_p,X_{p+1},\cdots ,X_{p+q}$ and $X_{p+q+r}$ needs to be found. Forecasting methods generally include single-step predictions, multistep predictions, and rolling predictions. Among them, the rolling prediction can reach a certain training sample and reflects the relationship between the output sample and predicted data. In fact, the rolling prediction starts from the single-step prediction, then feeds the output value back to the import as part of the input for predicting future values (see

Table 3. Process of rolling prediction.

Step	Input	Output (Prediction Result)
1	$X_p,X_{p+1},\cdots ,X_{p+q}$	$X_{p+q+1}$
2	$X_{p+1},X_{p+2},\cdots ,X_{p+q+1}$	$X_{p+q+2}$
⁝	⁝	⁝
$r$	$X_{p+q-1},X_{p+q},\cdots ,X_{p+q+r-1}$	$X_{p+q+r}$

).

3.4. Weight Determination

Hybrid forecasting is to assign several kinds of single prediction methods to different weights to form a comprehensive forecasting model. It can accurately and reasonably use the valuable information of a single forecasting model, better adapt to future changes, and reduce forecasting risk.

The key to hybrid forecasting is to properly determine the weight of various forecasting methods, and reasonable weights will improve the prediction accuracy greatly. Common weight determination methods include the arithmetic average method, variance reciprocal method, mean square reciprocal method, simple weighting method, and linear programming method. Among them, the arithmetic averaging method is suitable for treating each individual model equally if it is not known to each model. Usually, it is not optimal, and the sum of squared errors is large. The variance reciprocal method, mean square reciprocal method, and simple weighting method all have higher prediction accuracy than the arithmetic average method, but it is necessary to have a certain understanding of the prediction requirements in advance. All three methods above have one thing in common, that is, the variance is used to calculate the weight. The variance is the degree of fluctuation of the response variable above and below the mean. When the prediction result is not ideal and the error fluctuation is not large, the weight will still be large, which is unreasonable.

Therefore, this paper adopts the linear programming method, which determines the optimal weighting coefficient by taking the minimum absolute value of the combined prediction value as the objective function.

4. Case Study and Discussion

4.1. Data Collection

This paper collected the operating data of the Huilongwan site of China Auto Rental Inc. in Chongqing for 4 consecutive weeks (28 days). Each set of data includes the date of pick up, vehicle brand, car rental price, and pick-up location. According to the previous method, the model is divided into three grades according to the rental price and is set as a “presale lead time interval” every week. The survey data is sorted out to obtain three “presale lead time intervals” for each price grade model. The observable order quantity q and the presale opening status of each grade model are shown in

Table 4. Car rental order data.

Variable	${\mathbf{\Delta }}_1$							${\mathbf{\Delta }}_2$							${\mathbf{\Delta }}_3$
	1	2	3	4	5	6	7	1	2	3	4	5	6	7	1	2	3	4	5	6	7
$IB(1,j,t)$	11	5	7	13	3	9	4	6	3	5	3	9	12	5	10	7	14	10	14	12	14
$I(1,j,t)$	1	0	1	1	0	1	0	1	0	1	0	1	1	1	1	1	0	0	0	1	1
$IB(2,j,t)$	5	9	8	4	9	4	6	3	3	7	5	1	8	2	11	9	9	8	7	8	8
$I(2,j,t)$	0	1	1	0	1	0	1	1	0	1	1	0	1	1	1	0	1	0	1	1	0
$IB(3,j,t)$	3	1	4	2	3	1	1	3	1	2	3	5	3	2	4	2	2	1	2	2	3
$I(3,j,t)$	1	0	1	1	1	0	1	1	0	1	1	1	1	1	1	0	0	1	1	0	1

.

4.2. Applying IMSM

According to the rental price, the car type is classified as 1 under $40 per day, 2 within $40–70 per day, and 3 above $70 per day. Due to cost and time constraints, a random sampling method was used to conduct a questionnaire survey at some rental sites of China Auto Rental Inc. in Chongqing (see Figure 5).

The sample size of this RP survey was 400, and 363 questionnaires were collected, of which 320 were valid and 220 of which included buyup or buydown. Discretize the survey data, and let “very satisfied = 10”, “satisfied = 8”, “general = 5”, “not very satisfied = 3”, and “not satisfied = 0”. The influencing factor score, utilities, and possibilities corresponding to various choice behaviors can be obtained (

Table 5. Influencing factor scores, utilities, and possibilities.

Car Type Choice	Condition Attribute						Utility	Possibility
	${\mathit{c}}_3$	${\mathit{c}}_5$	${\mathit{c}}_6$	${\mathit{c}}_7$	${\mathit{c}}_8$	${\mathit{c}}_9$
1-0	7.14	6.64	6.74	6.33	6.42	6.71	6.74604	0.6156
1-2	6.35	6.18	5.22	6.35	5.58	5.26	6.00614	0.2938
1-3	4.34	4.97	4.86	5.45	5.42	4.61	4.82926	0.0905
2-0	5.64	6.35	6.33	5.48	5.77	6.13	5.95663	0.3547
2-1	5.89	6.87	6.24	5.87	6.15	5.29	6.14929	0.4301
2-3	5.11	6.05	6.11	4.93	5.25	5.11	5.45660	0.2152
3-0	6.92	6.85	6.83	6.32	6.54	6.77	6.76680	0.5558
3-1	5.89	4.46	5.66	5.89	6.12	6.03	5.48990	0.1550
3-2	6.36	6.15	6.29	5.36	6.02	6.15	6.11340	0.2892

).

Based on a week’s data from the International Convention & Exhibition Center, the biggest site of China Auto Rental Inc. in the Nan‘an district, EM, PD, MSM, and IMSM were used to repair the car rental demand, and the unconstrained repair results were compared with the real needs. The convergence conditions were $\theta =0.1$ and $\tau =0.6$, set in EM and PD (see

Table 6. Repair results and errors analysis of 1 week.

Variable	Estimation Method	1	2	3	4	5	6	7
$IB(1,j,t)$		45	46	53	51	43	49	47
$TU(1,j,t)$		49	46	53	51	47	49	49
$I(1,j,t)$		0	1	1	1	0	1	0
	EM	51.76 (5.63%)				51.59 (9.77%)		52.16 (6.45%)
	PD	50.11 (2.27%)				49.84 (6.04%)		50.7 (3.47%)
	MSM	49.82 (1.67%)				49.76 (5.87%)		50.04 (2.12%)
	IMSM	49.56 (1.14%)				47.68 (1.45%)		50.04 (2.12%)
$IB(2,j,t)$		24	24	24	23	20	20	26
$TU(2,j,t)$		24	24	24	23	22	20	26
$I(2,j,t)$		1	1	1	1	0	1	1
	EM					24.38 (10.82%)
	PD					23.5 (6.82%)
	MSM					23.27 (5.77%)
	IMSM					20.5 (6.82%)
$IB(3,j,t)$		12	14	12	12	10	15	15
$TU(3,j,t)$		13	14	12	12	13	15	15
$I(3,j,t)$		0	1	1	1	0	1	1
	EM	14.88 (14.46%)				14.64 (12.62%)
	PD	14 (10.82%)				13.63 (14.45%)
	MSM	12.85 (1.15%)				13.56 (4.31%)
	IMSM	13.12 (0.92%)				12.23 (5.92%)

Numbers in brackets are the errors of the estimated value relative to the real demands.

).

According to the results listed in the above table, when demand repair is not carried out, the average gap between constrained demand and actual demand is large, which is 10.1%. After various demand estimation methods are adopted, the average error of the obtained demand data is small, among which the average error of IMSM is the smallest, 3.06%.

4.3. Demand Estimation

Combined with the survey data, the abovementioned various algorithms are used to estimate the unconstrained demand of multiple car types. For EM and PD, let $\theta =0.1$ and $\tau =0.6$. In the PD algorithm, the difference of parameter $\tau$ affects the result of demand estimation and maximum iteration. The smaller $\tau$ is, the larger the iteration time is, and the higher the estimation of unconstrained demand becomes. When $\tau \le 0.35$, the estimation result of PD starts to be larger than EM. $\tau =0.6$ is selected here. For MSM and IMSM, the specific values of the parameters $\mu$ and $\sigma$ are listed in

Table 7. Estimation results comparison of EM, PD, MSM, and IMSM.

Method	Variable	${\mathbf{\Delta }}_1$							${\mathbf{\Delta }}_2$							${\mathbf{\Delta }}_3$
		1	2	3	4	5	6	7	1	2	3	4	5	6	7	1	2	3	4	5	6	7
	$IB(1,j,t)$	11	5	7	13	3	9	4	6	3	5	3	9	12	5	10	7	14	10	14	12	14
EM	$IU(1,j,t)$	11	12	7	13	12	9	12	6	10	5	10	9	12	5	10	7	18	15	18	12	14
PD	$IU(1,j,t)$	11	10	7	13	10	9	11	6	8	5	8	9	12	5	10	7	16	13	16	12	14
MSM		$\mu =10,\hspace{0.17em}\hspace{0.17em}\sigma =2.85$							$\mu =7.4,\hspace{0.17em}\hspace{0.17em}\sigma =3.05$							$\mu =10.75,\hspace{0.17em}\hspace{0.17em}\sigma =2.99$
	$SU(1,j,t)$	0	5	0	0	7	0	6	0	5	0	5	0	0	0	0	0	1	2	1	0	0
	$U(1,j,t)$	0	5	0	0	7	0	6	0	5	0	5	0	0	0	0	0	1	2	1	0	0
	$IU(1,j,t)$	11	10	7	13	10	9	10	6	8	5	8	9	12	5	10	7	15	12	15	12	14
IMSM	$CP(1,j,t)$	0	5	0	0	7	0	6	0	5	0	5	0	0	0	0	0	1	2	1	0	0
	$U(1,j,t)$	0	5	0	0	7	0	6	0	4	0	5	0	0	0	0	0	1	1	1	0	0
	$IU(1,j,t)$	11	10	7	13	10	9	10	6	7	5	8	9	12	5	10	7	15	11	15	12	14
	$IB(2,j,t)$	5	9	8	4	9	4	6	3	3	7	5	1	8	2	11	9	9	8	7	8	8
EM	$IU(2,j,t)$	11	9	8	11	9	11	6	3	8	7	5	7	8	2	11	12	9	11	7	8	11
PD	$IU(2,j,t)$	9	9	8	8	9	8	6	3	6	7	5	6	8	2	11	11	9	10	7	8	10
		$\mu =8,\hspace{0.17em}\hspace{0.17em}\sigma =1.41$							$\mu =5,\hspace{0.17em}\hspace{0.17em}\sigma =2.55$							$\mu =8.75,\hspace{0.17em}\hspace{0.17em}\sigma =1.71$
MSM	$SU(2,j,t)$	3	5	0	4	7	4	6	0	6	0	5	5	0	0	0	1	1	2	1	0	1
	$U(2,j,t)$	3	0	0	4	0	4	0	0	2	0	0	5	0	0	0	1	0	1	0	0	1
	$IU(2,j,t)$	8	9	8	8	9	8	6	3	5	7	5	6	8	2	11	10	9	9	7	8	9
IMSM	$CP(2,j,t)$	3	0	0	4	0	4	0	0	3	0	0	5	0	0	0	1	0	2	0	2	0
	$U(2,j,t)$	3	0	0	4	0	4	0	0	1	0	0	5	0	0	0	1	0	1	0	1	0
	$IU(2,j,t)$	8	9	8	8	9	8	6	3	4	7	5	6	8	2	11	10	9	9	7	9	7
	$IB(3,j,t)$	3	1	4	2	3	1	1	3	1	2	3	5	3	2	4	2	2	1	2	2	3
EM	$IU(3,j,t)$	3	4	4	2	3	4	1	3	5	2	3	5	3	2	4	4	4	1	2	4	3
PD	$IU(3,j,t)$	3	3	4	2	3	3	1	3	4	2	3	5	3	2	4	4	4	1	2	4	3
		$\mu =2.6,\hspace{0.17em}\hspace{0.17em}\sigma =1.14$							$\mu =3,\hspace{0.17em}\hspace{0.17em}\sigma =1.1$							$\mu =2.5,\hspace{0.17em}\hspace{0.17em}\sigma =1.29$
MSM	$SU(3,j,t)$	3	7	0	4	7	6	6	0	8	0	5	4	0	0	0	2	1	2	1	1	1
	$U(3,j,t)$	0	2	0	0	0	2	0	0	2	0	0	0	0	0	0	1	1	0	0	1	0
	$IU(3,j,t)$	3	3	4	2	3	3	1	3	3	2	3	5	3	2	4	3	3	1	2	3	3
IMSM	$CP(3,j,t)$	0	2	0	0	0	2	0	0	3	0	0	0	0	0	0	1	1	0	0	1	0
	$U(3,j,t)$	0	2	0	0	0	1	0	0	2	0	0	0	0	0	0	1	1	0	0	1	0
	$IU(3,j,t)$	3	3	4	2	3	2	1	3	3	2	3	5	2	3	4	3	3	1	2	3	3

for different vehicle types and cycles.

4.4. Demand Forecast

4.4.1. Holt–Winters Model Forecast

When making predictions, the first step is to calibrate the model. Due to the large amount of data, only the fitting process of car type 1 will be described in detail. Considering the characteristics of the Holt–Winters model, the initial value of cubic exponential smoothing is selected from the first value $IB(1,j,1)$. Therefore, these three sets of data are eliminated, and the remaining 28 sets of data are fitted. According to the minimum sum of squared residuals, the coefficients of $\alpha$, $\beta$, and $\gamma$ are determined to be 0.3, 0.9, and 0.1, respectively. The two curves obtained coincide with each other, which indicates that the exponential smoothing effect is good. The forecasting model is generally considered to be available when the error is within 10% (see Figure 6). The demands of different car types can be acquired. The specific prediction results are shown in

Table 8. Holt–Winters forecast results based on different estimation methods.

Method	Car Type	${\mathbf{\Delta }}_4$							Subtotal	Total
		1	2	3	4	5	6	7
EM	1	7	12	10	14	11	11	14	79	175
	2	7	8	8	7	8	18	15	71
	3	3	6	1	3	5	5	2	25
PD	1	11	3	10	10	8	13	12	67	168
	2	16	2	3	12	12	9	15	69
	3	1	6	2	3	6	5	9	32
MSM	1	12	2	5	7	11	13	13	63	155
	2	11	9	9	7	10	8	11	65
	3	5	3	5	3	2	4	5	27
IMSM	1	11	5	8	9	10	14	3	60	133
	2	10	9	9	8	12	5	1	54
	3	4	1	2	1	3	6	2	19

.

4.4.2. BP Neural Network Forecast

The BP neural network used in this paper is divided into three layers: the first layer is the input layer, with input variables of $i$, $j$, and $t$. The second layer is a hidden layer. In this paper, the same sample is used to train the models of different neuron numbers. According to the test error, it is reasonable to determine 17 hidden neurons. The hidden layer activation function uses the sigmoid function. The third layer is the output layer, with the demand as the output element, a total of one neuron, and the activation function uses a linear activation function.

After the BP neural network model is established, samples need to be trained. In general, the training sample needs to meet four characteristics: (1) there is a close functional relationship between the input and object variable, and the latter will change significantly due to the change of the former; (2) input variables are independent of each other, and other components cannot be accurately calculated using the components of the input variables; (3) the data to be predicted has a certain commonality with the sample data; (4) the sample capacity should have a certain scale, so that all the samples can reflect the mapping relationship between output variables and target variables.

In order to train and test the neural network model more reasonably, the unconstrained historical data repaired by EM is divided into 42 groups, among which seven groups of prediction samples and 35 groups of training samples are randomly selected by scrolling to select data. The training samples are used to train the neural network model, while the prediction samples are only used to test the feasibility and reliability of the forecasting model. The specific grouping situation is shown in

Table 9. Prediction samples and training samples.

No.	Car Type	Input Sample							Out Sample
1	1	13.00	11.49	9.00	11.49	6.00	9.74	5.00	9.74
2	1	12.00	5.00	10.00	7.00	17.13	15.04	17.13	12.00
3	2	9.00	8.00	10.31	9.00	10.30	6.00	3.00	7.38
4	2	3.00	7.38	7.00	5.00	6.96	8.00	2.00	11.00
5	2	8.00	2.00	11.00	11.06	9.00	10.61	7.00	8.00
6	2	2.00	11.00	11.06	9.00	10.61	7.00	8.00	10.61
7	3	3.00	4.26	2.00	3.00	5.00	3.00	2.00	4.00
8	1	11.00	11.57	7.00	13.00	11.49	9.00	11.49	6.00
9	1	11.57	7.00	13.00	11.49	9.00	11.49	6.00	9.74
10	1	7.00	13.00	11.49	9.00	11.49	6.00	9.74	5.00
11	1	11.49	9.00	11.49	6.00	9.74	5.00	9.74	9.00
12	1	9.00	11.49	6.00	9.74	5.00	9.74	9.00	12.00
13	1	11.49	6.00	9.74	5.00	9.74	9.00	12.00	5.00
14	1	6.00	9.74	5.00	9.74	9.00	12.00	5.00	10.00
15	1	9.74	5.00	9.74	9.00	12.00	5.00	10.00	7.00
16	1	5.00	9.74	9.00	12.00	5.00	10.00	7.00	17.13
17	1	9.74	9.00	12.00	5.00	10.00	7.00	17.13	15.04
18	1	9.00	12.00	5.00	10.00	7.00	17.13	15.04	17.13
19	1	5.00	10.00	7.00	17.13	15.04	17.13	12.00	14.00
20	2	10.38	9.00	8.00	10.31	9.00	10.30	6.00	3.00
21	2	8.00	10.31	9.00	10.30	6.00	3.00	7.38	7.00
22	2	10.31	9.00	10.30	6.00	3.00	7.38	7.00	5.00
23	2	9.00	10.30	6.00	3.00	7.38	7.00	5.00	6.96
24	2	10.30	6.00	3.00	7.38	7.00	5.00	6.96	8.00
25	2	6.00	3.00	7.38	7.00	5.00	6.96	8.00	2.00
26	2	7.38	7.00	5.00	6.96	8.00	2.00	11.00	11.06
27	2	7.00	5.00	6.96	8.00	2.00	11.00	11.06	9.00
28	2	5.00	6.96	8.00	2.00	11.00	11.06	9.00	10.61
29	2	6.96	8.00	2.00	11.00	11.06	9.00	10.61	7.00
30	3	3.00	3.48	4.00	2.00	3.00	3.48	1.00	3.00
31	3	3.48	4.00	2.00	3.00	3.48	1.00	3.00	4.26
32	3	4.00	2.00	3.00	3.48	1.00	3.00	4.26	2.00
33	3	2.00	3.00	3.48	1.00	3.00	4.26	2.00	3.00
34	3	3.00	3.48	1.00	3.00	4.26	2.00	3.00	5.00
35	3	3.48	1.00	3.00	4.26	2.00	3.00	5.00	3.00
36	3	1.00	3.00	4.26	2.00	3.00	5.00	3.00	2.00
37	3	4.26	2.00	3.00	5.00	3.00	2.00	4.00	3.89
38	3	2.00	3.00	5.00	3.00	2.00	4.00	3.89	3.89
39	3	3.00	5.00	3.00	2.00	4.00	3.89	3.89	1.00
40	3	5.00	3.00	2.00	4.00	3.89	3.89	1.00	2.00
41	3	3.00	2.00	4.00	3.89	3.89	1.00	2.00	3.89
42	3	2.00	4.00	3.89	3.89	1.00	2.00	3.89	3.00

.

Using prediction samples and training samples, the calculated values of forecasting models are checked separately. The neural network forecasting model has a high degree of fitness for the training samples, and the average error is 0.066 (see Figure 7). The model is reasonable and can be applied to later predictions. Demands of different car types can be obtained, which are shown in

Table 10. BP neural network forecast results based on different estimation methods.

Method	Car Type	${\mathbf{\Delta }}_4$							Subtotal	Total
		1	2	3	4	5	6	7
EM	1	12	14	10	9	8	16	16	85	189
	2	12	8	10	9	10	13	13	75
	3	5	3	4	3	2	5	7	29
PD	1	13	13	10	10	8	13	17	84	184
	2	11	12	8	9	7	14	12	73
	3	5	3	2	3	3	5	6	27
MSM	1	12	13	12	9	8	13	15	82	174
	2	11	11	6	9	5	13	13	68
	3	5	3	2	2	3	4	5	24
IMSM	1	11	12	8	9	7	11	16	74	163
	2	10	12	6	8	5	12	14	67
	3	4	4	2	1	3	3	5	22

.

4.4.3. Hybrid Model

The combination weight coefficients are selected by a linear programming method. When the weight is 0.7 for the Holt–Winters model and 0.3 for the BP neural network, the sum of the squared errors is the smallest. Using the hybrid forecasting model, the demands of different car types can be obtained, as shown in

Table 11. Hybrid forecasting results based on different estimation methods.

Method	Car Type	${\mathbf{\Delta }}_4$							Subtotal	Total
		1	2	3	4	5	6	7
EM	1	10	13	10	11	9	14	15	82	181
	2	10	11	9	8	8	15	12	73
	3	4	4	2	3	3	5	5	26
PD	1	12	9	10	10	8	13	15	77	176
	2	13	8	6	10	9	12	13	71
	3	3	4	2	3	4	5	7	28
MSM	1	12	9	9	8	9	13	14	74	165
	2	11	10	7	8	7	11	12	66
	3	5	3	3	2	3	4	5	25
IMSM	1	11	10	8	9	8	12	12	2	154
	2	10	11	7	8	7	10	10	63
	3	4	3	2	1	3	4	4	21

. Hybrid forecasting can neutralize both models, making the predictions more accurate, increasing the smoothness, and forecasting accuracy, so it is presented below.

4.5. Discussion

As described previously, the demand data was estimated by four unconstrained estimation methods, and the three forecast models were used to predict the demand in the fourth presale lead time interval; then, 12 sets of demand data were obtained. In contrast to the real demand in the fourth presale lead time interval (see

Table 12. Real demand in the fourth presale lead time interval.

Car Type								Subtotal	Total
	1	2	3	4	5	6	7
1	11	8	11	9	8	10	12	69	147
2	9	7	7	8	7	9	11	58
3	3	2	2	2	3	4	4	20

), the relative errors of the various results can be calculated (see

Table 13. Relative errors of forecast demand for each car type and all car types.

Forecast Model	Car Type	Unconstrained Estimation Method
		EM	PD	MSM	IMSM
Holt–Winters	1	12.31%	−3.62%	−10.10%	−12.10
	2	20.69%	17.24%	8.62%	7.47%
	3	6.00%	47.45%	32.50%	−6.67%
BP Neural Network	1	23.19%	21.74%	18.84%	7.27%
	2	29.31%	25.86%	17.24%	15.52%
	3	36.00%	35.00%	20.00%	10.00%
Hybrid Model	1	18.84%	11.59%	7.25%	1.45%
	2	25.86%	22.41%	13.79%	8.66%
	3	24.00%	40.00%	25.00%	5.00%
Holt–Winters BP Neural Network Hybrid Model	All types	14.96%	11.56%	3.06%	9.50%
		28.57%	25.17%	18.37%	10.89%
		23.12%	19.72%	12.24%	4.76%

).

From the real demand and the prediction results (as shown in Table 12 and Table 13), the following can be seen.

For different car types, from the perspective of the predicted effect, comparing the Holt–Winters model with the BP neural network model, the former has a greater volatility, while the latter is more stable. From the perspective of the adaptability of the unconstrained estimation model, when the unconstrained estimation data obtained by MSM is used, the prediction results of the three methods are relatively stable, and the hybrid forecasting model can obtain the best results. The relative errors of the three models are respectively: 1.45%, 8.66%, and 5.00% (see data in the box of Table 13).

For the forecasted total, the Holt–Winters model is better than the BP neural network model in predicting the effect. From the perspective of the adaptability of the unconstrained estimation model, the Holt–Winters model has the best prediction effect when using MSM to repair the data, while the three methods can all obtain better prediction results when using IMSM, among which hybrid model has the smallest error.

Considered comprehensively, when using MSM for unconstraining repair data, the relative error of the Holt–Winters model for predicting total demand is 3.06%. However, looking at the corresponding subtype prediction results, we find that the relative error of the third type is 32.5%. When IMSM is used for unconstrained repair, the relative error of the total demand using the hybrid model is 4.76%. The relative errors of the predicted results for the three models are 1.45%, 8.66%, and 5.00%. Therefore, by using the IMSM to estimate historical constraint requirements, the hybrid forecasting is optimal for future demand prediction.

5. Conclusions

To fully use the information of historical car rental data, a two-stage joint approach was proposed for predicting multitype car rental demand. The method considers the customer’s choice behavior to improve the spill model, which can estimate the unconstrained demand effectively and uses the hybrid forecasting model to predict future short-term demand.

For historical car rental demand data, the repaired data can more accurately reflect customer demand. In order to prove the effectiveness of the method, EM, PD, and MSM were compared in the unconstrained demand estimation stage, and the Holt–Winters, BP neural network, and hybrid models were used for comparison in the future demand prediction stage. For different car types, different models can form different prediction results in two stages. The comprehensive calculation results show that the proposed method is superior to other combinations in terms of stability and effectiveness of prediction.

Specifically, using the case study, for the repair of the demand data of each car type, IMSM has obvious advantages compared with the other three methods, and the average error of the example test is 3.06%. For the prediction of future demand of each car type, based on the calculation results of different unconstrained estimation models, the performances of the Holt–Winters model, the BP neural network model, and the hybrid model are quite different. The hybrid forecasting model has the best effect, and the relative error of the predicted results is within the acceptable range for different car types or overall needs. Therefore, the case study shows that the proposed method outperforms other methods. It can be regarded as a new way to help car rental companies to predict customer demand. There are many factors affecting customers’ car rental behavior, and with the application of mobile apps, the way of car reservation is changing. In the future, more influencing factors, such as electric vehicles, urban traffic restriction policies, and road traffic conditions, should be taken into account. In addition, the different distribution characteristics of car rental demand and the combined optimization of forecasting methods will also be interesting research topics.