1 Introduction

In queueing frameworks, that system is considered as efficient which limits the long-run cost. Now a days due to the developing anxiety of clients for not waiting in queues during the briefest span, estimation of minimum waiting time is also done by few authors to develop a robust system. However in optimization, concept of convexity plays a vital role since it can guarantee the presence of global optimal solution. A local minima of a convex function on a convex feasible region is a global minima. In many papers related to optimization of cost or waiting time of queueing models, the optimal solution is obtained using optimization tools but the sufficient optimality conditions are not verified. In consequence we have tried to design a model which closely relates to the system developed for facilitating mankind and optimized not only the cost of the system but awaiting time as well. Also, we have shown the optimization problems as a convex programming problem by discussing the convexity of the objective functions using Hessian matrix test and tried to obtain the global minima of the minimization problem.

Queueing models with retrials find their application in most modern computer networks that include Wide Area Network (WAN) protocol, Frame Relay, Local Area Network (LAN) protocol, TCP/IP and X.25 as well. Mohammadi et al. (2014) have designed a reliable healthcare network where the service provider is erratic with queues of finite capacity. Gao et al. (2020) inspected a queue with retrials along with distinct failures in the model. Also if the server is unoccupied and faces breakdown (passive breakdown), it gets a delayed repair due to unidentified fault whereas if the server is occupied and faces breakdown (active breakdown) then the repair process starts immediately. They have given a promising application of the model in packet-switching network. Before sending forward, the messages are sloted into IP packets which are then sent from a source to destination through the router. Here source host and the destination host are arrivals and servers, respectively. More applications of retrial queueing models with breakdown can be seen in Upadhyaya (2014), Kim and Kim (2016), Taleb and Aissani (2016), Lan and Tang (2020), Lee et al. (2020), Dragieva and Phung-Duc (2020).

1.1 Literature review

Study of stochastic models is done by many researchers. Choudhury and Tadj (2009) extended the classical M/G/1 model with secondary service, breakdowns and delayed repair where erratic server randomly breaks-down while serving the clients and further evaluated the busy and awaiting time along with reliability indices. On the same lines, Taleb and Aissani (2016) also incorporated the impatient behaviour of the clients with secondary services while Saggou et al. (2017) focused mainly on recurrent clients.

Now we discuss a well known characteristic of clients called balking (resisting) where the arriving batch may get annoyed and resist to join the queue if the worker is occupied or on holiday. A queue having bulk arrival with retrials, balking together with altered vacations was studied by Upadhyaya (2014). In addition, the upshot of some basic parameters on long run probabilities and mean orbit size is also studied. Seeing the recent developments, Upadhyaya and Kushwaha (2020) elaborated \(M^X/G/1\) queueing model with retrials incorporating the impatience of clients with their feedback and delayed repair. They considered alternate vacation scheme and validated the analytical results using a soft computing technique known as ANFIS (Adaptive Neuro-Fuzzy Interface System) and then built a cost effective system.

Multi-objective optimization is of interest of researchers working in this field. In multiobjective optimization more than one objective function is to be optimized simultaneously so that the Pareto optimal solution can be obtained for the problem. Verma (1986) first introduced the concept of multiobjective optimization in queueing systems. Kahag et al. (2019) did a bi-objective optimization of M/M/m model for a hub allocation problem in traffic systems using a multi objective invasion weed optimization (MOIWO). Most recently, Wu and Yang (2021) considered a two phase heterogeneous service model. They first did a single objective optimization using Canonical PSO and then formulated a bi-objective cost optimization model for the system and awaiting time simultaneously. More studies on multi-objective optimization in queueing models can be read in the work of Mohammadi et al. (2014), Tavakkoli-Moghaddam et al. (2017) and Pourmohammadi et al. (2021).

Although lot of work is done in obtaining the optimal cost of the system for queueing models but only a handful of authors have discussed the concept of convexity. Zhang et al. (1997) studied M/G/1 continuous model with two distinct vacations and developed the average cost function with optimal threshold policies. Later, Zhang (2006) evaluated the cost function’s convexity for the model studied in Zhang et al. (1997). Sherman et al. (2009) have examined an erratic queue with retrial where normal and orbit queue have infinite capacity and the clients join the orbit if the server fails. They allow both active as well as idle breakdowns in their system. They have obtained the optimal repair and retrial rates and have done the convexity analysis in detail.

1.2 Focus of our study

Motivated by the pertinence of such models, M/G/1 queue with retrials under active and passive breakdown services, delayed repair for passive faults with balking is studied in this work. The main focus is to attain the global minima of the cost function and further getting the global Pareto optimal solution for the bi-objective problem. Application of the developed model can be seen in health care systems, supermarkets, call centres, mobile and computer network organisations. This model is closely related with the smooth functioning of the system during some internal faults (virus attack, electricity failures etc.). Our model firmly reflects the vaccination system developed where effective regular and delayed repair can assume an essential part in conveying quality help rapidly under decreased expense. The manufactures and decision-makers can efficiently apply the outcomes of this model in developing system management policies. In addition to its functional pertinence, the model presented interesting numerical properties that permits the examination by their own doing.

1.3 Layout of the paper

The remaining article is composed as follows: Sect. 2 includes the basic notations and preliminaries that are to be considered throughout. Section 3 depicts queueing model with a practical application. The steady state study of the model is referred in Sect. 4 where the governing and boundary equations, along with the probability generating function (pgf) of the system and performance indices are evaluated. Section 5 elaborates the effect of some important parameters on the queue length with ANFIS validation. The single and multi-objective optimization of the model is discussed in Sect. 6 where in sub-segment 6.1, the expected cost function is fabricated to accomplish the ideal estimations of certain parameters and get the optimal cost value using PSO. Later in sub-segment 6.2, the bi-objective optimization where two objectives: the awaiting time in the system and the cost function are optimized and lastly Sect. 7 incorporate concluding remarks along with future scope of the paper.

2 Preliminaries and notations

2.1 Preliminaries

For \({\bar{\omega }},{\bar{\Omega }}\in R^n\), following notations for inequalities will be followed as:

$$\begin{aligned} \begin{array}{ll} {\bar{\omega }}\leqq {\bar{\Omega }}\iff {\bar{\omega }}_j\leqq {\bar{\Omega }}_j,&{}\quad j=1,\cdots ,n~;\\ {\bar{\omega }}\le {\bar{\Omega }}\iff {\bar{\omega }}\leqq {\bar{\Omega }},&{}\quad {\bar{\omega }}\ne {\bar{\Omega }};\\ {\bar{\omega }}< {\bar{\Omega }}\iff {\bar{\omega }}_j< {\bar{\Omega }}_j,&{}\quad j=1,\cdots ,n. \end{array} \end{aligned}$$

Definition 1

Any function \(T: R^n\rightarrow R\) is convex if for \({\tilde{x}}\in R^n\) we have \(T[k {\tilde{x}}+(1-k){\tilde{x}}]\leqq kT({\tilde{x}})+(1-k)T({\tilde{x}})\) true \(\forall \) k where \(0\leqq k \leqq 1\). The function T is strictly convex if the above inequality holds strictly for \(x\ne {\tilde{x}},\) and \(0<k <1\).

Definition 2

If all partial derivatives of T exists and are continuous over the domain, then the Hessian matrix M is a square \(n\times n\) matrix arranged as

$$\begin{aligned} M= \left[ \begin{array}{llll} \frac{{\partial ^2}{T}}{\partial {x_1}^2}&{} \frac{{\partial ^2}{T}}{\partial x_1,\partial x_2} &{} \cdots &{} \frac{{\partial ^2}{T}}{\partial x_1,\partial x_n}\\ \frac{{\partial ^2}{T}}{\partial x_2,\partial x_1}&{} \frac{{\partial ^2}{T}}{\partial x_2} &{} \cdots &{} \frac{{\partial ^2}{T}}{\partial x_2,\partial x_n}\\ \vdots &{} \vdots &{} \ddots &{} \vdots \\ \frac{{\partial ^2}{T}}{\partial x_n,\partial x_1}&{} \frac{{\partial ^2}{T}}{\partial x_n,\partial x_2} &{} \cdots &{} \frac{{\partial ^2}{T}}{\partial x_n}\\ \end{array} \right] ~. \end{aligned}$$

Assume S to be a nonempty open convex set in \(R^n\), and let \(T: R^n\rightarrow R\) be two times differentiable on S. Then T is said to be convex iff the matrix M is positive semidefinite at each point in S. For more details refer Bazaraa et al. (2005).

Consider an optimization problem

where \(x\in X\) is an open subset of \(R^n\) and \(T:X\rightarrow R^k\) are differentiable on X. If \(k=1\), problem is called a scalar otherwise multi objective optimization problem.

Definition 3

(Global Optimal point) If in MP \(k=1\), any point \({\tilde{x}}\in R^n\) is an optimal point if \(~\exists \) no \(x\in X\) s.t. \(T(x)< T({\tilde{x}})\).

Definition 4

(Global Pareto Optimality) A point \({\tilde{x}}\) is a global efficient/Pareto optimal solution of MP if \(\exists \) no \(x\in X\) s.t. \(T(x)\le T({\tilde{x}}).\)

Obviously any global Pareto solution is local Pareto optimal but converse is not always true. The converse holds for a convex multi-objective programming problem. Also, the multi-objective problem of optimization is convex if all the objectives and feasible region is convex. For more details refer Miettinen (1998).

2.2 Basic notations

$$\begin{aligned} \begin{array}{ll} \lambda &{} :\text {Poisson arrival rate of the clients.} \\ \alpha &{} :\text {Exponential retrial rate of the clients.}\\ b (\bar{b}=1-b), &{} :\text {If the system is occupied then the client's joining (resisting) }\\ &{} \text {probability from the orbit.} \\ \delta , \theta &{} :\text {Exponential breakdowns, rate for passive and active}\\ &{} \text {breakdowns, respectively.} \\ B, R, G &{} :\text {Random variable denoting general distribution of the serving time, }\\ &{} \text {repair time of active and passive breakdowns, respectively.}\\ B(x_1), R(x_1), G(x_1) &{} :\text {Cumulative density function (cdf) of the serving time,}\\ &{} \text {repair time of active and passive breakdowns, respectively.}\\ b(x_1), r(x_1), g(x_1) &{} :\text {Probability generating function (pgf) of the serving time,}\\ &{} \text {repair time of active and passive breakdowns, respectively.}\\ \beta _{i}, v_{i}, \mu _{i} &{} :i^{th} \text {moment of serving time, repair time of active }\\ &{} \text {and passive breakdowns, respectively.} \\ \beta (x_1), v(x_1), \mu (x_1) &{} :\text {Conditional complementary rates for serving time,} \\ &{} \text {repair time of active and passive breakdowns, respectively.}\\ \end{array} \end{aligned}$$

3 Model representation and Practical illustration

In this segment, the mathematical model with stability condition is included along with a practical application of the model in real life scenario. In sub-segment. 3.1, a detailed model description is given with a schematic diagram as in Fig. 1 for the readers. sub-segment 3.2 provides a practical application using example for better understanding of the model. Further in sub-segment.  3.3 Markov chain is defined and the system’s condition to be stable is stated.

3.1 Model description

We are considering an erratic retrial line with two distinct active and passive breakdowns and delayed repair because of the latter. The basic assumptions for the model are:

  • Arrivals and balking: The arrival rate of a client is \(\lambda \) following Poisson process. If the event is that the server is occupied or is uncertain, then, at that point the arriving client may get annoyed and resist joining the queue with b probability or may exit from the system with \(\bar{b}\) probability.

  • Service and retrial policy: The entering client, on seeing the server unoccupied gets the service straight away. Else on seeing the server occupied or inoperative, the clients will have to retry later from a virtual orbit following FCFS principle. The client on the top of the virtual track (orbit) retries to get served with exponential rate \(\alpha \), if the service provider is unoccupied.

  • Breakdowns and repairs: The server may face active and passive breakdowns in occupied and unoccupied states, respectively. In the active breakdown, the server just instantly starts repairing. While in passive breakdown, the unnoticed server stays inactive until a client enters and then only the server enters delayed repair state. In both the cases the service of the client continues (starts) its service once the server is ready to serve.

Fig. 1
figure 1

Schematic diagram of the model

Full size image

We assume the conditional complementary rates when the service provider is occupied, for repair time of active and passive failures, respectively, as

$$\begin{aligned} \beta (y)=\frac{b(y)}{\bar{B}(y)}, v(y)=\frac{r(y)}{\bar{R}(y)}, \mu (y)=\frac{g(y)}{\bar{G}(y)}. \end{aligned}$$

For the remaining paper, if we consider CDF as F(y), where y is any variable, then its complement is \(\bar{F}(y)=1-F(y)\), its Laplace-Stieltjes transform (LST) is \({\tilde{F}}(\kappa )=\int _{0}^{\infty }e^{-{\kappa y}}dF(y)\) and the Laplace transform (LT) of \(\bar{F}(y)\) is \(\bar{F}^{*}(\kappa )=\int _{0}^{\infty }e^{-{\kappa y}}(1-F(y))dy\) which gives \(\bar{F}^{*}(\kappa )=\frac{1-{\tilde{F}}(y)}{\kappa }\).

3.2 Practical examples

A practical example for our retrial queue can be seen if a client goes to a bank to deposit cash or get the passbook updated. We assume that their is only one employ (service provider) in the bank who is assigned to do the job of cash deposit and passbook updates. The clients arriving the bank may get immediate service if the computer system (server) is idle. Else they will have to wait in a queue and try again after some time (retrial) for their service if the computer is already occupied with the earlier present client. Now the retrying clients may get restless and may want to exit from the bank without getting their work done (balk) if they don’t have much time to wait. In that case the client may leave the bank (system) without service. We know that the bank’s employ have all the data and information of all the clients in their computer system and being a machinery product it may face some faults (breakdown). On a new day initially the system is idle and if there is any fault in the computers at that point (passive breakdown) the bank’s employ won’t know it. As soon as the client enters the bank, the employ starts the computer to serve him and only then he will know about the fault in the computer and will send it to get repaired (delayed repair). The time interval between point when the passive fault occurred and the point at which the client entered is termed as delayed period. The client that arrived during delayed period will be provided service once the repair work is done. Another possibility is that while the employer is serving the client the computer system may get a virus attack and hence may breakout (active breakdown). In this situation, it will have to be send for immediate repair. This bank scenario can be easily modelled with our retrial model having two faults and delayed repair.

The model closely fits to the recent Covid-19 scenario also where the main focus of all the countries and governments is to vaccinate maximum people as early as possible. It is not an easy task for a country like India with population of around 140 crores. Thus the whole vaccination process is divided into two stages. First stage includes the booking of the vaccination slots. Individuals have to book their respective slots to get vaccinated through an app like Cowin app in India. Our model fits into this situation where the individuals who are supposed to be vaccinated can be treated as patients and the website/app/computer system is the service provider where the slots are to be booked. The system may face some faults in the idle state (passive breakdown) whose repair will start once the patient starts the system for booking. Thus, there will be a delay in the repair for a passive breakdown (delayed repair). When the system is repaired and is ready to work upon, the patient will continuously try to book the vaccination slot. In India, the booking slots were restricted on daily basis. As lakhs of people try booking at the same time, the server could be busy at particular instances and so the patients have to retry for the booking after an interval of time (retrials). Some of the patients might get impatient and may want to resist from retrying, so they may leave the system without booking the slot (balking). Sometimes due to continuous retrials and busy server, the system may get heated up or may face a virus attack and may breakdown in the middle of the booking process (active breakdown). In this situation the system will be repaired immediately and then the service will be accomplished. Finally, when a slot would be booked for an individual, the service is completed. The second phase includes the vaccination center allotted while slot booking. Here the servers are the nurses assigned for vaccination purpose. The queue of patients wait outside the center for their turn. We aim to bring down the cost of the system along with the awaiting time of the patients. All the components of our model fits in this real-life situation and thus we try to get the optimal repair and retrial rates so that one can book the slot at the earliest and get vaccinated with minimum system cost and awaiting time.

3.3 Stability condition

We consider \(S_G\) to be the generalized time interval to serve the clients from starting to the end of the service, where CDF is \(S_G(y)\) and LST is \({\tilde{S}}_G(y)\). If we consider the possibility of active faults in the service, we get \({\tilde{S}}_G(y)={\tilde{B}}(y+\theta (1-{\tilde{R}}(y)))\), which gives \(E[S_{G}]=\beta _{1}(1+\theta v_{1})=\beta _{1}^{*}\) and \(E[S_{G}^{2}]=\beta _{2}(1+\theta v_{1})^{2}+\theta \beta _{1} v_{2}=\beta _{2}^{*}\).

Next we define

$$\begin{aligned} \begin{aligned} a_k&=\int _{0}^{\infty }\frac{(\lambda bt)^k}{k!}e^{-\lambda bt}dG(t), ~~~k=0,1,2,\ldots ,\\ h_k&=\int _{0}^{\infty }\frac{(\lambda bt)^k}{k!}e^{-\lambda bt}dS_G(t), ~~~k=0,1,2,\ldots , \end{aligned} \end{aligned}$$

where \(a_k\) denotes the likelihood of k clients entering the orbit if the service provider is under repair of passive fault, \(h_k\) denotes the likelihood of k clients joining the virtual track during generalized service time. Defining \(c_k=\sum _{i=0}^{k}a_i h_{k-i}, ~k\geqq 0\) where \(c_k\) denotes the likelihood of k clients entering in orbit when the service provider is in generalized service and under the passive repair time.

Let \(\Delta (x_3)=\lambda b(1-x_3)\). Then we define

$$\begin{aligned} \begin{aligned} H(x_3)&=\sum _{k=0}^{\infty }{x_3}^k h_k ={\tilde{B}}(\Delta (x_3)+\theta (1-{\tilde{R}}(\Delta (x_3)))),\\ H'(1)&=\frac{\mathrm{d}{H(x_3)}}{\mathrm{d}t}\Big |_{{x_3}=1}^{}=\lambda b\beta _{1}^{*}=\rho ,\\ H''(1)&=\frac{{\mathrm{d}^2}{H(x_3)}}{{\mathrm{d}t}}\Big |_{{x_3}=1}^{}=(\lambda b)^{2}\beta _{2}^{*},\\ A(x_3)&=\sum _{k=0}^{\infty }{x_3}^k a_k ={\tilde{G}}(\Delta (x_3)),\\ A'(1)&=\frac{\mathrm{d}{A(x_3)}}{{\mathrm{d}t}}\Big |_{{x_3}=1}^{}=\lambda b\mu _{1}=\rho _{1},\\ A''(1)&=\frac{{\mathrm{d}^2}{A(x_3)}}{{\mathrm{d}t}}\Big |_{{x_3}=1}^{}=(\lambda b)^{2}\mu _{2},\\ C(x_3)&=\sum _{k=0}^{\infty }{x_3}^k c_k =A(x_3)H(x_3). \end{aligned} \end{aligned}$$

Considering \(T_k\) (\(T_0 =0\)) as the time epoch when the kth client exits from the system, \(O_k=O(T_k)\) be the size of the virtual track when the Kth exit takes place, thus \(O_{k} for k>=0\) is a Markov chain process with state space given by \({\mathbb {N}}\). The theorem stated below gives the necessary and sufficient condition for the stable system.

Theorem 3.1

The Markov chain \(\{O_k, K\geqq 0\}\) is stable iff \(\rho +\frac{\delta \rho _1}{\lambda b+\alpha +\delta }<\frac{\alpha }{\lambda b+\alpha }\).

Proof

One can see the proof of Theorem 3.1 from Gao et al. (2020) and can establish the prove for the above theorem on same lines. \(\square \)

4 Scrutinizing the steady state

We use supplementary variable technique (SVT) to procure the steady state solution along with the pgf of system size in this sub-plot. At a particular time t, the Markov process defines the system’s state as \(\{N'(t), S(t), \tau _{1}(t), \tau _{2}(t), \tau _{4}(t), t\geqq 0\}\), where \(N'(t)\) is the patient’s number and S(t) shows the server’s state which is as follows:

$$\begin{aligned} S(t)= \left\{ \begin{array}{ll} 0, &{} \text {unoccupied server,} \\ 1, &{} \text {occupied server,} \\ 2, &{} \text {repair of server's active fault,} \\ 3, &{} \text {server's delayed repair period,} \\ 4, &{} \text {repair of passive breakdown of the server,} \\ \end{array} \right. \end{aligned}$$

where when \(S(t)=1\), \(\tau _{1}(t)\) gives the passed service time; \(S(t)=2\), \(\tau _{2}(t)\) gives the passed repair time of active breakdown; \(S(t)=4\), \(\tau _{4}(t)\) gives the passed repair time of passive breakdown;

The steady state probabilities and their densities are defined as follows:

$$\begin{aligned} \begin{aligned} Q_{n,i}&=P(N'=n, ~S=i)\\&=\lim _{t\rightarrow \infty } Q_{n,i}(t), ~~~n\geqq 0, ~~~i=0,3,\\ Q_{n,i}(x_1)dx_1&=P(N'=n, ~S=i, ~~x_1<\tau _{i}\leqq x_1+dx_1)\\&=\lim _{t\rightarrow \infty } Q_{n,i}(t,x_1)dx_1, ~~~n\geqq 0, ~~~i=1,4, ~~~x_1\geqq 0,\\ Q_{n,2}(x_1,x_2)dx_1dx_2&=P(N'=n, ~S=2, ~~x_1<\tau _{1}\leqq x_1+dx_1, ~~x_2<\tau _{4}\leqq x_2+dx_2)\\&=\lim _{t\rightarrow \infty } Q_{n,i}(t,x_1,x_2)dx_1dx_2, ~~~n\geqq 0, ~~~x_1,x_2\geqq 0,\\ \end{aligned} \end{aligned}$$

4.1 The governing and boundary equations

Using the SVT we proceed with the basic equations defining steady states as:

$$\begin{aligned}&(\lambda +\delta )Q_{0,0}=\int _{0}^{\infty } Q_{0,1}(x_1)\beta (x_1)dx_1, \end{aligned}$$
(1)
$$\begin{aligned}&(\lambda +\delta +\alpha )Q_{n,0}=\int _{0}^{\infty } Q_{n,1}(x_1)\beta (x_1)dx_1, ~~~n\geqq 1, \end{aligned}$$
(2)
$$\begin{aligned}&\frac{d}{dx_1} Q_{n,1}(x_1)=-(\lambda b+\theta +\beta (x_1)) Q_{n,1}(x_1)+\lambda b Q_{n-1,1}(x_1)\nonumber \\&\quad +\int _{0}^{\infty } Q_{n,2}(x_1,x_2)v(x_2)dx_2, ~~~ n\geqq 0, ~~~x_1>0, \end{aligned}$$
(3)
$$\begin{aligned}&\frac{d}{dx_1} Q_{n,2}(x_1,x_2)=-(\lambda b+v(x_2)) Q_{n,2}(x_1,x_2)\nonumber \\&\quad +\lambda b Q_{n-1,2}(x_1,x_2), ~~ n\geqq 0, ~~x_1,x_2>0, \end{aligned}$$
(4)
$$\begin{aligned}&\lambda b Q_{0,3}=\delta Q_{0,0}, \end{aligned}$$
(5)
$$\begin{aligned}&(\lambda b+\alpha ) Q_{n,3}=\delta Q_{n,0}, ~~~n\geqq 1, \end{aligned}$$
(6)
$$\begin{aligned}&\frac{d}{dx_1} Q_{n,4}(x_1)=-(\lambda b+\mu (x_1)) Q_{n,4}(x_1)+\lambda b Q_{n-1,4}(x_1), ~~ n\geqq 0, ~~x_1>0 . \end{aligned}$$
(7)

The boundary conditions that are to be used in solving Eqs. (1)–(7) are:

$$\begin{aligned}&Q_{n,1}(0)=\lambda b Q_{n,0}+\alpha Q_{n+1,0} +\int _{0}^{\infty } Q_{n,4}(x_1)\mu (x_1)dx_1, ~~~ n\geqq 0, \end{aligned}$$
(8)
$$\begin{aligned}&Q_{n,2}(x_1,0)=\theta Q_{n,1}, ~~ n\geqq 0, ~~x_1>0, \end{aligned}$$
(9)
$$\begin{aligned}&Q_{n,4}(0)=\lambda b Q_{n,3}+\alpha Q_{n+1,3}, ~~ n\geqq 0, \end{aligned}$$
(10)

where the normalizing condition is

$$\begin{aligned} \begin{aligned} \sum _{n=0}^{\infty }(Q_{n,0}+Q_{n,3})&+\sum _{n=0}^{\infty }\int _{0}^{\infty }(Q_{n,1}(x_1)+Q_{n,4}(x_1))dx_1\\&+\sum _{n=0}^{\infty }\int _{0}^{\infty }\int _{0}^{\infty }Q_{n,2}(x_1,x_2)dx_1dx_2 =1.\\ \end{aligned} \end{aligned}$$
(11)

4.2 Probability generating function for the system states

We next define the generating functions that are to be used during further calculations as:

$$\begin{aligned} \begin{aligned} Q_{i}(x_3)&=\sum _{n=0}^{\infty }Q_{k,i}{x_3}^{n}, ~~~i=0,3,\\ Q_{i}(x_1,x_3)&=\sum _{n=0}^{\infty }Q_{n,i}(x_1){x_3}^{n}, ~~~i=1,4,\\ Q_{2}(x_1,x_2,x_3)&=\sum _{n=0}^{\infty }Q_{n,2}(x_1,x_2){x_3}^{n}.\\ \end{aligned} \end{aligned}$$

Theorem 4.1

The generating functions for the joint stationary distribution of the states of the service provider are as follows:

$$\begin{aligned} \begin{aligned} Q_{0}(x_3)&=\frac{\alpha (\lambda b+\alpha )(x_3-H(x_3))-\delta \alpha C(x_3)(1-x_3)}{\bar{D}(x_3)} Q_{0,0} , \\ Q_{1}(x_1,x_3)&=\frac{N_{a}(x_3)}{\bar{D}(x_3)}(1-B(x_1))\alpha Q_{0,0} \times exp\{-(\Delta (x_3)+\theta (1-{\tilde{R}}(\Delta (x_3))))x_1\}, \\ Q_{2}(x_1,x_2,x_3)&=\theta Q_{1}(x_1,x_2)(1-R(x_2)) exp\{-\Delta (x_3)x_2\}, \\ Q_{3}(x_3)&=\frac{\delta Q_{0,0}}{\lambda b+\alpha }\times \left[ \frac{\alpha }{\lambda b}+\frac{\alpha (\lambda b+\alpha )(x_3-H(x_3))-\delta \alpha C(x_3)(1-x_3)}{\bar{D}(x_3)}\right] , \\ Q_{4}(x_1,x_3)&=\frac{N_{b}(x_3)}{\bar{D}(x_3)}(1-G(x_1))exp\{-\Delta (x_3)x_1\}\delta \alpha Q_{0,0} , \\ \end{aligned} \end{aligned}$$

where

$$\begin{aligned} \begin{aligned} Q_{0,0}&=\frac{\lambda b}{\alpha }\frac{(\alpha +\lambda b)(\alpha +\lambda b+\delta )}{(\delta +\lambda b)(\alpha +\lambda b+\delta )+\delta (\alpha +\lambda b)\rho _{1}}\left( \frac{\alpha }{\alpha +\lambda b}-\rho -\frac{\delta \rho _{1}}{\alpha +\lambda b+\delta } \right) , \\ N_{a}(x_3)&=\delta A(x_3) (\alpha +\lambda bx_3+(\alpha +\lambda b+\delta )(x_3-1))-(\alpha +\lambda b)(\delta +\Delta (x_3)),\\ N_{b}(x_3)&=(\alpha +\lambda bx_3)(1-H(x_3))+(\alpha +\lambda b+\delta )(x_3-1),\\ \bar{D}(x_3)&=(\alpha +\lambda b)[(\alpha +\lambda b+\delta )x_3-(\alpha +\lambda bx_3)H(x_3)]-\delta (\alpha +\lambda bx_3)C(x_3).\\ \end{aligned} \end{aligned}$$

Proof

Using Eqs. (1)–(10) and doing some algebraic manipulations, the proof of the theorem can be obtained easily. \(\square \)

Theorem 4.2

If there is a service provider which is occupied, under repair of active or passive fault, then the marginal pgfs of the size of orbit are as follows:

$$\begin{aligned} \begin{aligned} Q_{1}(x_3)&=\frac{N_{a}(x_3)}{\bar{D}(x_3)}\frac{1-H(x_3)}{\Delta (x_3)+\theta (1-{\tilde{R}}(\Delta (x_3)))}\alpha Q_{0,0},\\ Q_{2}(x_3)&=\frac{N_{a}(x_3)}{\bar{D}(x_3)}\frac{1-H(x_3)}{\Delta (x_3)+\theta (1-{\tilde{R}}(\Delta (x_3)))}\frac{1-{\tilde{R}}(\Delta (x_3))}{\Delta (x_3)}\theta \alpha Q_{0,0},\\ Q_{4}(x_3)&=\frac{N_{b}(x_3)}{\bar{D}(x_3)}\frac{1-A(x_3)}{\Delta (x_3)}\delta \alpha Q_{0,0}.\\ \end{aligned} \end{aligned}$$

Proof

We already have marginal probabilities for unoccupied and delayed repair state from the above theorem. So, we now need to get the marginal probabilities for the busy and both the repair states. Thus evaluate \(Q_{1}(x_3)=\int _{0}^{\infty }Q_{1}(x_1,x_3)dx_1\), \(Q_{2}(x_3)=\int _{0}^{\infty }\int _{0}^{\infty }Q_{2}(x_1,x_2,x_3)dx_1dx_2\) and \(Q_{4}(x_3)=\int _{0}^{\infty }Q_{4}(x_1,x_3)dx_1\) using Theorem 4.1 and hence the proof can be obtained. \(\square \)

Theorem 4.3

The PGFs \(\phi (x_3)\) and \(\psi (x_3)\) of the client’s number in the orbit and system respectively are

$$\begin{aligned} \begin{aligned} \phi (x_3)&=\frac{(x_3-1)(\alpha +\delta +\lambda b)(\delta +\lambda b)+\delta (\alpha +\lambda b)[C(x_3)-H(x_3)]}{\bar{D}(x_3)}\frac{\alpha }{\lambda b}Q_{0,0}~~,\\ \psi (x_3)&=\frac{\delta C(x_3) [(\delta +\lambda b)(x_3-1)+(\alpha +\lambda b)x_3]-(\alpha +\lambda b)(\delta +\Delta (x_3))H(x_3)}{\bar{D}(x_3)}\frac{\alpha }{\lambda b}Q_{0,0}~.\\ \end{aligned} \end{aligned}$$

Proof

Let us suppose that \(\phi (x_3)=E[{x_3}^{N_{0}}]\) and \(\psi (x_3)=E[{x_3}^{N_{1}}]\) where \(N_{0}\) and \(N_{1}\) denotes the customer’s number in the orbit and system respectively. Then using \(\phi (x_3)=\sum _{n=0}^{\infty }Q_{j}(x_3)\) and \(\psi (x_3)=Q_{0}(x_3)+{x_3}Q_{1}(x_3)+{x_3}Q_{2}(x_3)+Q_{3}(x_3)+{x_3}Q_{4}(x_3)\), proof of the theorem can be obtained. \(\square \)

4.3 Performance indices

The aim of this sub-segment is to provide the important performance indices of the queueing system using the results of sub-segment 4.2.

Theorem 4.4

(A) Taking into account the steady state conditions,:

  • \(Q_{a}\)- probability of service provider being unoccupied is

    $$\begin{aligned} Q_{a}=Q_{0}(1)=\frac{\lambda b[\delta +(\alpha +\lambda b)(1-\rho )]}{(\delta +\lambda b)(\alpha +\lambda b+\delta )+\delta \rho _{1}(\alpha +\lambda b)}. \end{aligned}$$

  • \(Q_{b}\)- probability that the service provider is occupied is

    $$\begin{aligned} Q_{b}=Q_{1}(1)=\frac{\lambda b(\alpha +\lambda b)+\delta ((\alpha +\lambda b)\rho _{1}+\delta +\alpha +2\lambda b)}{(\delta +\lambda b)(\alpha +\lambda b+\delta )+\delta \rho _{1}(\alpha +\lambda b)}\lambda b \beta _{1}. \end{aligned}$$

  • \(Q_{c}\)- probability of service provider being under repair of active fault is

    $$\begin{aligned} Q_{c}=Q_{2}(1)=\theta v_{1}Q_{b}. \end{aligned}$$

  • The probability \(Q_{d}\) that the service provider is under delayed repair is

    $$\begin{aligned} Q_{d}=Q_{3}(1)=\frac{(\alpha +\lambda b+\delta )(1-\rho ) -\lambda b\rho -\delta \rho _{1}}{(\delta +\lambda b)(\alpha +\lambda b+\delta )+\delta \rho _{1}(\alpha +\lambda b)}\delta . \end{aligned}$$

  • The probability \(Q_{e}\) that the service provider is under repair of passive fault is

    $$\begin{aligned} Q_{e}=Q_{4}(1)=\frac{\delta +(\alpha +\lambda b)(1-\rho )}{(\delta +\lambda b)(\alpha +\lambda b+\delta )+\delta \rho _{1}(\alpha +\lambda b)}\delta \rho _{1}. \end{aligned}$$

(B) The mean orbit (\(L_{0}\)) and system size (\(L_{1}\)) are

$$\begin{aligned} \begin{aligned} L_{0}&=\frac{N^{''}(1)}{2N^{'}(1)}-\frac{\bar{D}^{''}(1)}{2\bar{D}^{'}(1)}=Q_f+Q_g,\\ L_{1}&=Q_b+Q_c+Q_e+L_{0}. \end{aligned} \end{aligned}$$
(12)

if

$$\begin{aligned} \begin{aligned} N&=(x_3-1)(\alpha +\delta +\lambda b)(\delta +\lambda b)+\delta (\alpha +\lambda b)(C(x_3)-H(x_3)),\\ N^{'}(1)&=(\alpha +\delta +\lambda b)(\delta +\lambda b)+\delta (\alpha +\lambda b)\rho _{1},\\ N^{''}(1)&=2\delta \lambda b\rho _{1}+\delta (\alpha +\lambda b)((\lambda b)^{2}\mu _{2}+2\rho \rho _{1}),\\ \bar{D}^{'}(1)&=(\alpha +\delta +\lambda b)(\alpha +\lambda b)\left[ \frac{\alpha }{\alpha +\lambda b}-\rho -\frac{\delta \rho _{1}}{\alpha +\lambda b+\delta }\right] ,\\ \bar{D}^{''}(1)&=-(\alpha +\lambda b)\left[ 2\lambda b\rho +(\alpha +\lambda b)((\lambda b)^{2}\beta ^{*}_{2})\right] ,\\&\quad -\delta \left[ 2\lambda b(\rho +\rho _{1})+(\alpha +\lambda b)((\lambda b)^{2}\mu _{2}+2\rho \rho _{1}+(\lambda b)^{2}\beta ^{*}_{2})\right] ~.\\ \end{aligned} \end{aligned}$$

(C) The mean awaiting time in the orbit (\(W_{0}\)) and system (\(W_{1}\)) are

$$\begin{aligned} \begin{aligned} W_{0}&=\frac{L_{0}}{\lambda }~,\\ W_{1}&=\frac{L_{1}}{\lambda }~. \end{aligned} \end{aligned}$$
(13)

Proof

  1. (A)

    Substituting \({x_3}=1\) in Theorem 4.1 and 4.2 and directly calculating using L’Hopital’s differentiation rule, we can prove it.

  2. (B)

    Using \(L_{0}=E[N_{0}]= d{\phi (x_3)}{t}_{{x_3}=1}^{}\) and \(L_{1}=E[N_{1}]= d{\psi (x_3)}{t}_{{x_3}=1}^{}\) the mean orbit and system size can be obtained. (C) From Little’s theorem, the results can be established easily. \(\square \)

4.4 Special cases

Here we present some particular cases. It can be seen that under certain conditions, the system size of our model reduces to the expression obtained for previously studied models.

  1. (i)

    When \(b=1\) (No balking), we observe that the model reduces to the model studied by Gao et al. (2020) and the results coexist.

  2. (ii)

    When \(\delta =0\) (No passive breakdown), we observe the change as the model changes to continuous queue with constant retrial and active breakdown and the results coincide with that of Jin-ting (2006) on considering \(A(u)=1-e^{-\alpha u} ; u > 0\).

  3. (iii)

    When \(\delta =0, \theta =0, b=1\) (No breakdowns, no repair and no balking), the results of our model coincide with the model studied by Gomes Corral (1999).

  4. (iv)

    When \(b=1\), \(\theta =0\) (No balking and no repair), the results of our model coincide with those of Taleb and Aissani (2016) on considering no persistent, impatient clients and preventive maintenance in their model.

5 Numerical instance

For the numerical example, we assume that repair times and service times follow Erlangian distribution. The repair time for passive failure has mean \(\mu _{1}=5/\mu \) and variance \(\mu _{2}=5/\mu ^2\). The time to repair the active fault has mean \(v_{1}=2/v\) and variance \(v_{2}=2/v^2\). The service time has mean \(\beta _{1}=2/s\) and variance \(\beta _{2}=2/s^2\). The retrial time and passive and active breakdown values are given as: \(\lambda =0.15, ~b=0.25, ~\alpha =2, ~\delta =0.2, ~\theta =0.1, ~\mu = 5, ~v=5, ~s=0.75\).

A soft computing technique, ANFIS is a robust tool to discover significant outcomes that are applicable in everyday crowding situation. This technique helps in finding estimated solution for the measures whose definite outcomes are otherwise difficult to obtain. In ANFIS technique, we handle the parameters (\(\lambda , \alpha , s\) and \(\theta \)) as linguistic variables which are executed for 4 epochs each. The linguistic values for all the parameters are defined as follows: less, moderate, high and extreme. We have used Gaussian membership function to present all the linguistic variables as shown in Figs. 2b, 3b, 4b and 5b for \(\lambda , \alpha , s\) and \(\theta \), respectively. In Figs. 2a, 3a, 4a and 5a, the solid line shows the analytical results whereas dotted lines with solid symbol show the ANFIS results. It is clear from the figures that the analytical results for the model coincide with those of the neuro fuzzy results attained using ANFIS technique.

By applying this, we aim to present the effect on mean system size by changing \(\lambda , \alpha , s\) and \(\theta \). It can be observed that \(L_1\) increases exponentially with rising retrials and decreases in logarithmic manner with increasing service rate. With the increase in \(\lambda \) and \(\theta \), \(L_1\) rises gradually in almost same manner. Thus the system designers and decision makers should make efforts and choose these parameters in such a way that their respective system becomes robust.

Fig. 2
figure 2

a \(L_1\) vs \(\lambda \) for analytic and ANFIS results; b membership function for \(\lambda \)

Full size image
Fig. 3
figure 3

a \(L_1\) vs \(\alpha \) for analytic and ANFIS results; b membership function for \(\alpha \)

Full size image
Fig. 4
figure 4

a \(L_1\) vs s for analytic and ANFIS results; b membership function for s

Full size image
Fig. 5
figure 5

a \(L_1\) vs \(\theta \) for analytic and ANFIS results; b membership function for \(\theta \)

Full size image

6 Cost optimization

Practically, the total operating cost of a system assumes a vital part in the investigation of many industrial systems. System architects and supervisors are typically keen on limiting the working expense per unit time to make the system more productive and profitable. More work for optimization of the cost function can be seen in Gao et al. (2020), Lan and Tang (2020), Upadhyaya and Kushwaha (2020). In this segment, the optimal outline of the retrial queue including distinct failures and delayed repair is addressed. The analysis of a cycle for the system can be seen in detail in Gao et al. (2020) from where using the argument of an alternating renewal process, we get

$$\begin{aligned} E[\psi ]=\frac{1}{(\lambda +\delta )Q_{0,0}} \end{aligned}$$

where

$$\begin{aligned} \begin{aligned} Q_{0,0}=\frac{\lambda b}{\alpha }\frac{BC}{AB+\delta C\rho _1}\left( \frac{\alpha }{C}-\rho -\frac{\delta }{B}\rho _1\right) ,\\ A=\lambda b+\delta , ~~~~B= \alpha +\lambda b+\delta ~~~and ~~~~C= \alpha +\lambda b.\\ \end{aligned} \end{aligned}$$

To exhibit the relevance of the outcomes acquired in the past conversation, we foster an expected operating cost function for the queueing model per unit time which is defined as

$$\begin{aligned} \begin{aligned} TC&= C_hL_1 + C_uQ_a + C_bQ_b + C_aQ_c + C_dQ_d + C_pQ_e + C_s\frac{1}{E[\psi ]}, \\ TC&= C_hL_1 + C_uQ_a + C_bQ_b + C_aQ_c + C_dQ_d + C_pQ_e + C_s(\lambda +\delta )Q_{0,0}.\\ \end{aligned} \end{aligned}$$
(14)

where the cost symbols corresponding to the cost function (per unit time) are as listed below:

\(C_h\) : Cost for holding each client in the system,

\(C_u\) : Cost of service provider being unoccupied,

\(C_b\) : Cost of service provider being occupied,

\(C_a\) : Cost of service provider being under repair for an active breakdown,

\(C_d\) : Cost of service provider being under delayed period,

\(C_p\) : Cost of service provider being under repair for a passive breakdown,

\(C_s\) : Set up cost for the busy cycle that’s fixed.

6.1 Single objective optimization

Here we aim to minimize the system cost TC and try to find the optimal repair and retrial rates. Substituting the values of \(L_1, Q_a, Q_b, Q_c, Q_d, Q_e\) and \(Q_{0,0}\) in Eq. (14), we get

$$\begin{aligned} \begin{aligned} TC=&\,TC(\alpha , v)= C_h \left( \frac{2\delta \lambda b\rho _1+\delta C((\lambda b)^2\mu _2+2\rho \rho _1)}{2[AB+\delta C\rho _{1}]}\right. \\&\left. +\frac{C[2\lambda b\rho +C(\lambda b)^2{\beta _2}^*]+\delta [2\lambda b(\rho +\rho _1)+C((\lambda b)^2\mu _2+2\rho \rho _1+(\lambda b)^2{\beta _2}^*)]}{2BC[\frac{\alpha }{C}-\rho -\frac{\delta }{B}\rho _1]} \right. \\&\left. +\frac{\lambda bC+\delta (C\rho _1+\delta +\alpha +2\lambda b)}{AB+\delta C\rho _{1}}\lambda b\beta _1 +\theta v_1\frac{\lambda bC+\delta (C\rho _1+\delta +\alpha +2\lambda b)}{AB+\delta C\rho _{1}}\lambda b\beta _1 \right. \\&\left. +\delta \rho _{1}\frac{\delta +C(1-\rho )}{AB+\delta C\rho _{1}}\right) + C_u \frac{\lambda b[\delta +C(1-\rho )]}{AB+\delta C\rho _{1}}\\&+ C_b \frac{\lambda bC+\delta (C\rho _1+\delta +\alpha +2\lambda b)}{AB+\delta C\rho _{1}}\lambda b\beta _1+ C_a\theta v_1\frac{\lambda bC+\delta (C\rho _1+\delta +\alpha +2\lambda b)}{AB+\delta C\rho _{1}}\lambda b\beta _1\\&+ C_d \delta \frac{B(1-\rho )-\lambda b\rho -\delta \rho _1}{AB+\delta C\rho _{1}}+ C_p\delta \rho _{1}\frac{\delta +C(1-\rho )}{AB+\delta C\rho _{1}}\\&+ C_s(\lambda +\delta )\frac{\lambda b}{\alpha }\frac{BC}{AB+\delta C\rho _1}\left( \frac{\alpha }{C}-\rho -\frac{\delta }{B}\rho _1\right) .\\ \end{aligned} \end{aligned}$$
(15)

It can be clearly seen that Eq. (15) is a nonlinear function of decision variables \(\alpha \) and v that are continuous. Thus, the minimization problem to get the minimum system cost can be described as follows:

$$\begin{aligned} \begin{aligned} (CP) ~~~~~~~~~~~~~~~~~~~~TC(\alpha ^*, v^*)&= \min TC(\alpha , v) \end{aligned} \end{aligned}$$
(16)

where \(S:\{(\alpha , v): \alpha , v\in (0.2,10)\}\subset R^2\).

Clearly S is a convex feasible set. The certainty that the normal expense function is complex implies that we cannot use traditional slope based methodologies to minimize. Consequently, we utilize a heuristic calculation to manage the issue of advancement in Eq. (16). The Particle Swarm Optimization algorithm (PSO) presented by Kennedy amb Eberhart (1995) is an adaptable optimization strategy that is performed to solve nonlinear objective functions. This algorithm works on the principle that the best positioned solution attract other possible values to get the optimal result in a particular search space.

For the optimization purpose the basic assumptions for critical parameters are:

\(\lambda =0.15, ~b=0.25, ~s=0.75, ~\delta =0.2, ~\theta =0.1, ~\mu = 5\).

Furthermore we have considered four distinct cost sets for the evaluation of minimum cost which are given as:

Cost Set-I:

Cost Set-II:

Cost Set-III:

Cost Set-IV:

Under such conditions, we study the changing pattern of the cost function for changing values of \(\lambda \) and \(\theta \) in Table 1. Clearly the rise in the cost of the system with a rise in these parameters can be observed. Using the MATLAB code for PSO algorithm, the total optimal cost of the system obtained for cost set-I is \({\mathbf{TC}(\alpha , \mathbf{v})=\$ \mathbf{17}}\) for \({(\alpha ^*, \mathbf{v}^*)={\textbf {(0.9439, 9.9995)}}}\). The convergence of the objective function using PSO is shown in Fig. 6.

Table 1 Using PSO, the effect of \(\lambda \) and \(\theta \) is obtained on the system cost
Full size table
Fig. 6
figure 6

Convergence curve using PSO

Full size image

Now we discuss the concept of convexity of the cost function. Here using the Hessian matrix test we have shown that the objective function proposed in the problem (CP) is convex and hence is a convex programming problem (CPP). Thus the optimal cost obtained above is a global minimum. Considering Eq. (14) as:

$$\begin{aligned} \begin{aligned} TC(\alpha , v)&= C_h(Q_b+Q_c+Q_e+Q_f+Q_g) + C_uQ_a + C_bQ_b \\&\quad + C_aQ_c+ C_dQ_d+ C_pQ_e + C_s(\lambda +\delta )Q_{0,0}.\\ \end{aligned} \end{aligned}$$
(17)

The second order partial derivatives of the cost function Eq. (17) are

$$\begin{aligned} \frac{{\partial ^2}{TC}}{{\partial \alpha }^2}= & {} C_h\left( \frac{{\partial ^2}{Q_b}}{{\partial \alpha }^2}+\frac{{\partial ^2}{Q_c}}{{\partial \alpha }^2}+\frac{{\partial ^2}{Q_e}}{{\partial \alpha }^2}+\frac{{\partial ^2}{Q_f}}{{\partial \alpha }^2}+\frac{{\partial ^2}{Q_g}}{{\partial \alpha }^2}\right) + C_u\frac{{\partial ^2}{Q_a}}{{\partial \alpha }^2} \nonumber \\&\quad + C_b\frac{{\partial ^2}{Q_b}}{{\partial \alpha }^2}+ C_a\frac{{\partial ^2}{Q_c}}{{\partial \alpha }^2}+ C_d\frac{{\partial ^2}{Q_d}}{{\partial \alpha }^2}+ C_p\frac{{\partial ^2}{Q_e}}{{\partial \alpha }^2} + C_s(\lambda +\delta )\frac{{\partial ^2}{Q_{0,0}}}{{\partial \alpha }^2}, \nonumber \\ \frac{{\partial ^2}{TC}}{\partial \alpha \partial v}= & {} C_h\left( \frac{{\partial ^2}{Q_b}}{\partial \alpha \partial v}+\frac{{\partial ^2}{Q_c}}{\partial \alpha \partial v}+\frac{{\partial ^2}{Q_e}}{\partial \alpha \partial v}+\frac{{\partial ^2}{Q_f}}{\partial \alpha \partial v}+\frac{{\partial ^2}{Q_g}}{\partial \alpha \partial v}\right) + C_u\frac{{\partial ^2}{Q_a}}{\partial \alpha \partial v} \nonumber \\&\quad + C_b\frac{{\partial ^2}{Q_b}}{\partial \alpha \partial v}+ C_a\frac{{\partial ^2}{Q_c}}{\partial \alpha \partial v}+ C_d\frac{{\partial ^2}{Q_d}}{\partial \alpha \partial v}+ C_p\frac{{\partial ^2}{Q_e}}{\partial \alpha \partial v} + C_s(\lambda +\delta )\frac{{\partial ^2}{Q_{0,0}}}{\partial \alpha \partial v}, \nonumber \\ \frac{{\partial ^2}{TC}}{{\partial v}^2}= & {} C_h\left( \frac{{\partial ^2}{Q_b}}{{\partial v}^2}+\frac{{\partial ^2}{Q_c}}{{\partial v}^2}+\frac{{\partial ^2}{Q_e}}{{\partial v}^2}+\frac{{\partial ^2}{Q_f}}{{\partial v}^2}+\frac{{\partial ^2}{Q_g}}{{\partial v}^2}\right) + C_u\frac{{\partial ^2}{Q_a}}{{\partial v}^2} \nonumber \\&\quad + C_b\frac{{\partial ^2}{Q_b}}{{\partial v}^2}+ C_a\frac{{\partial ^2}{Q_c}}{{\partial v}^2}+ C_d\frac{{\partial ^2}{Q_d}}{{\partial v}^2}+ C_p\frac{{\partial ^2}{Q_e}}{{\partial v}^2} + C_s(\lambda +\delta )\frac{{\partial ^2}{Q_{0,0}}}{{\partial v}^2}, \end{aligned}$$
(18)

where

$$\begin{aligned}&\frac{{\partial ^2}{Q_a}}{{\partial \alpha }^2} =\frac{{2(\lambda b + \delta (1 + {\rho _1}))[\lambda b\delta \rho A]}}{[AB+\delta C\rho _{1}]^3}; \quad \frac{{\partial ^2}{Q_b}}{{\partial \alpha }^2} =\frac{{\partial ^2}{Q_c}}{{\partial \alpha }^2}=0;\\&\frac{{\partial ^2}{Q_d}}{{\partial \alpha }^2} =\frac{{2\delta (\lambda b + \delta (1 + {\rho _1}))[{\delta ^2}(\rho + {\rho _1}) + \lambda b\delta \rho ]}}{[AB+\delta C\rho _{1}]^3}; ~~~~~~\frac{{\partial ^2}{Q_e}}{{\partial \alpha }^2} ={\rho _1}\frac{{\partial ^2}{Q_d}}{{\partial \alpha }^2};\\&\frac{{\partial ^2}{Q_f}}{{\partial \alpha }^2} =\frac{[\lambda b+\delta (1+\rho _1)][(\lambda b)^2\mu _2+2\rho \rho _1]\delta A+2\lambda b\delta \rho _1[\lambda b+\delta (1+\rho _1)]^2}{[AB+\delta C\rho _{1}]^3}; \\&\frac{{{\partial ^2}{Q_g}}}{{\partial {\alpha ^2}}} = \frac{\begin{aligned} {\left[ {\alpha B - \rho BC - \delta C{\rho _1}} \right] ^2}{(\lambda b)^2}\beta _2^* - (1 - \rho )\left[ {\alpha B - \rho BC - \delta C{\rho _1}} \right] \\ \times \left\{ {C\left[ {2\lambda b\rho + C{{(\lambda b)}^2}\beta _2^*} \right] + \delta \left[ {2\lambda b(\rho + {\rho _1}) + C({{(\lambda b)}^2}{\mu _2} + 2\rho {\rho _1} + {{(\lambda b)}^2}\beta _2^*)} \right] } \right\} \\ - \left[ {\alpha B - \rho BC - \delta C{\rho _1}} \right] \left[ {\alpha + B - 2\rho C - \delta (\rho + {\rho _1})} \right] \left\{ \begin{aligned} 2\lambda b\rho + 2C{(\lambda b)^2}\beta _2^* \\ + \delta \left( {{{(\lambda b)}^2}{\mu _2} + 2\rho {\rho _1} + {{(\lambda b)}^2}\beta _2^*)} \right) \\ \end{aligned} \right\} \\ + {\left[ {\alpha + B - 2\rho C - \delta (\rho + {\rho _1})} \right] ^2}\left\{ \begin{aligned} C\left[ {2\lambda b\rho + C{{(\lambda b)}^2}\beta _2^*} \right] \\ + \delta \left[ {2\lambda b(\rho + {\rho _1}) + C({{(\lambda b)}^2}{\mu _2} + 2\rho {\rho _1} + {{(\lambda b)}^2}\beta _2^*)} \right] \\ \end{aligned} \right\} \\ \end{aligned} }{{{{\left\{ {BC\left( {\frac{\alpha }{C} - \rho - \frac{{\delta {\rho _1}}}{B}} \right) } \right\} }^3}}}; \\&\frac{{\partial ^2}{Q_{0,0}}}{{\partial \alpha }^2} = \lambda b(\lambda + \delta )\left\{ \left( {\frac{{BC}}{{\alpha (AB + \delta C{\rho _1})}}} \right) \left( {\frac{{-2\lambda b}}{{{C^3}}} + \frac{-2\delta {\rho _1}}{{B^3}}} \right) + \frac{Y\left( \frac{\alpha }{C} - \rho - \frac{\delta \rho _1}{B} \right) }{(\alpha ^4)(AB + \delta C\rho _1)^4}\right. \\&\quad \left. + 2\left( {\frac{{\alpha (AB + \delta C{\rho _1})(2C + \delta ) - BCX}}{{{\alpha ^2}{{(AB + \delta C{\rho _1})}^2}}}} \right) \left( {\frac{{\lambda b}}{{{C^2}}} + \frac{{\delta {\rho _1}}}{{{B^2}}}} \right) \right\} ;\\ \text {such that}\\&X=[A(\alpha +B)+\delta (\alpha +C)\rho _1]\\&Y={\alpha ^2}{{(AB + \delta C{\rho _1})}^2}[(AB + \delta C{\rho _1})(2\lambda b + 4\alpha + \delta ) + (2\lambda b + 2\alpha + \delta )X(\alpha - 1) - BC(2A + 2\delta {\rho _1})] \\&-[\alpha (AB + \delta C{\rho _1})(2C + \delta ) - BCX] \times [2\alpha {(AB + \delta C{\rho _1})^2} + 2{\alpha ^2}(AB + \delta C{\rho _1})X]\\ {\vspace{6pt}}&\frac{{{\partial ^2}{Q_a}}}{{\partial \alpha \partial v}} = \frac{{2{{(\lambda b)}^2}{\beta _1}\theta \left[ {\left( {AB + \delta {\rho _1}C} \right) - C(\lambda b + \delta (1 + \rho ))} \right] }}{{{v^2}{{\left( {AB + \delta {\rho _1}C} \right) }^2}}};~~~~\frac{{{\partial ^2}{Q_b}}}{{\partial \alpha \partial v}} = 0; \\&\frac{{{\partial ^2}{Q_c}}}{{\partial \alpha \partial v}} = \frac{{ - 2\lambda b{\beta _1}\theta (\lambda b + \delta (1 + \rho ))\left[ {\left( {AB + \delta {\rho _1}C} \right) - (\lambda bC + \delta (C{\rho _1} + \lambda b + B))} \right] }}{{{v^2}{{\left( {AB + \delta {\rho _1}C} \right) }^2}}}; \\&\frac{{{\partial ^2}{Q_d}}}{{\partial \alpha \partial v}} = \frac{{2\theta \lambda b\delta {\beta _1}\left( {AB + \delta {\rho _1}C} \right) - \delta (\lambda b + \delta (1 + \rho ))\left[ {\lambda b{\beta _1}\left( {\frac{{2\theta }}{{{v^2}}}} \right) (B + \lambda b)} \right] }}{{{v^2}{{\left( {AB + \delta {\rho _1}C} \right) }^2}}}; \\&\frac{{{\partial ^2}{Q_e}}}{{\partial \alpha \partial v}} = \frac{{2\lambda b\theta {\beta _1}\delta {\rho _1}\left[ {\left( {AB + \delta {\rho _1}C} \right) - C(\lambda b + \delta (1 + \rho ))} \right] }}{{{v^2}{{\left( {AB + \delta {\rho _1}C} \right) }^2}}}; \\&\frac{{{\partial ^2}{Q_{0,0}}}}{{\partial \alpha \partial v}} = (\lambda + \delta )2\theta {(\lambda b)^2}{\beta _1}\left\{ {\frac{{\alpha \left( {AB + \delta {\rho _1}C} \right) (B + C) - BC[A(\alpha + B) + \delta {\rho _1}(\alpha + C)]}}{{{v^2}{\alpha ^2}{{\left( {AB + \delta {\rho _1}C} \right) }^2}}}} \right\} ; \\&\frac{{{\partial ^2}{Q_f}}}{{\partial \alpha \partial v}} = \frac{{ - 2\lambda b\theta {\beta _1}\delta {\rho _1}\left[ {\left( {AB + \delta {\rho _1}C} \right) - C(\lambda b + \delta (1 + \rho ))} \right] }}{{{v^2}{{\left( {AB + \delta {\rho _1}C} \right) }^2}}}; \\&\frac{{{\partial ^2}{Q_g}}}{{\partial \alpha \partial v}} = \frac{\begin{aligned} \left( {{{{4\lambda b\theta {\beta _1}BC}} {{{v^2}}}}} \right) \left\{ \begin{aligned} 2[(\alpha + B) - \rho (B + C) - \delta {\rho _1}]N' + 2BC\left( {\frac{\alpha }{C} - \rho - \frac{{\delta {\rho _1}}}{B}} \right) \frac{{\partial N'}}{{\partial \alpha }} \\ - \left( {{{{4\lambda b\theta {\beta _1}(B + C)}} {{{v^2}}}}} \right) N - \frac{{\partial N}}{{\partial \alpha }}\left( {{{{4\lambda b\theta {\beta _1}BC}} {{{v^2}}}}} \right) \\ \end{aligned} \right\} \\ - 2[(\alpha + B) - \rho (B + C) - \delta {\rho _1}]\left\{ {2BC\left( {\frac{\alpha }{C} - \rho - \frac{{\delta {\rho _1}}}{B}} \right) N' - \left( {{{{4\lambda b\theta {\beta _1}BC}} {{{v^2}}}}} \right) N} \right\} \\ \end{aligned} }{{{{\left\{ {2BC\left( {\frac{\alpha }{C} - \rho - \frac{{\delta {\rho _1}}}{B}} \right) } \right\} }^3}}} \end{aligned}$$

such that

$$\begin{aligned} \begin{aligned}&N = C\left[ {2{{(\lambda b)}^2}{\beta _1}\left( {1 + \frac{{2\theta }}{v}} \right) + C{{(\lambda b)}^2}\left( {{\beta _2}{{\left( {1 + \frac{{2\theta }}{v}} \right) }^2} + \theta {\beta _1}\frac{2}{{{v^2}}}} \right) } \right] \\&+ \delta \left[ {2{{(\lambda b)}^2}{\beta _1}\left( {1 + \frac{{2\theta }}{v}} \right) + 2\lambda b{\rho _1} + C\left( \begin{aligned} {(\lambda b)^2}{\mu _2} + 2{(\lambda b)^2}{\beta _1}\left( {1 + \frac{{2\theta }}{v}} \right) {\rho _1} \\ + {(\lambda b)^2}\left( {{\beta _2}{{\left( {1 + \frac{{2\theta }}{v}} \right) }^2} + \theta {\beta _1}\frac{2}{{{v^2}}}} \right) \\ \end{aligned} \right) } \right] , \\&N' = {(\lambda b)^2}BC\left[ {4\theta {\beta _2}\left( {1 + \frac{{2\theta }}{v}} \right) \left( {\frac{{ - 1}}{{{v^2}}}} \right) - \frac{{4\theta {\beta _1}}}{{{v^3}}}} \right] + 4\lambda b\theta {\beta _1}\left( {\frac{{ - 1}}{{{v^2}}}} \right) \left( {\lambda bB + \delta {\rho _1}C} \right) , \\&\frac{{\partial N}}{{\partial \alpha }} = 2\lambda b\rho + 2C{(\lambda b)^2}\beta _2^* + \delta ({(\lambda b)^2}{\mu _2} + 2\rho {\rho _1} + {(\lambda b)^2}\beta _2^*), \\&\frac{{\partial N'}}{{\partial \alpha }} = {(\lambda b)^2}\left[ {4\theta {\beta _2}\left( {1 + \frac{{2\theta }}{v}} \right) \left( {\frac{{ - 1}}{{{v^2}}}} \right) - \frac{{4\theta {\beta _1}}}{{{v^3}}}} \right] (B + C) + 4\lambda b\theta {\beta _1}\left( {\frac{{ - 1}}{{{v^2}}}} \right) \left( {\lambda b + \delta {\rho _1}} \right) . \\&\frac{{{\partial ^2}{Q_a}}}{{\partial {v^2}}} = \frac{{ - 4{{(\lambda b)}^2}{\beta _1}\theta C}}{{{v^3}{\left( {AB + \delta {\rho _1}C} \right) }}};~~~~~\frac{{{\partial ^2}{Q_b}}}{{\partial {v^2}}} = 0; \\ {\vspace{8pt}}&\frac{{{\partial ^2}{Q_c}}}{{\partial {v^2}}} = \frac{{4\lambda b{\beta _1}\theta (\lambda bC + \delta (C{\rho _1} + \lambda b + B))}}{{{v^3}\left( {AB + \delta {\rho _1}C} \right) }}; \\&\frac{{{\partial ^2}{Q_d}}}{{\partial {v^2}}} = \frac{{ - 4\lambda b\theta {\beta _1}\delta (\lambda b + B)}}{{v^3}{\left( {AB + \delta {\rho _1}C} \right) }};~~~~~\frac{{{\partial ^2}{Q_e}}}{{\partial {v^2}}} = \frac{{ - 4\lambda b\delta \theta {\beta _1}C{\rho _1}}}{{{v^3}\left( {AB + \delta {\rho _1}C} \right) }}; \\&\frac{{{\partial ^2}{Q_{0,0}}}}{{\partial {v^2}}} = - 4\theta \lambda b{\beta _1}\frac{{(\lambda + \delta )}}{{{v^3}}}\frac{{\lambda b}}{\alpha }\frac{{BC}}{{\left( {AB + \delta {\rho _1}C} \right) }};~~~~~\frac{{{\partial ^2}{Q_f}}}{{\partial {v^2}}} = \frac{{4\lambda b\theta {\beta _1}\delta {\rho _1}C}}{{{v^3}\left( {AB + \delta {\rho _1}C} \right) }}; \\&\frac{{{\partial ^2}{Q_g}}}{{\partial {v^2}}} = \frac{{2BC\left( {\frac{\alpha }{C} - \rho - \frac{{\delta {\rho _1}}}{B}} \right) {N_1}^{\prime \prime } - {N_1}4BC\theta \lambda b{\beta _1}\left( {\frac{{ - 2}}{{{v^3}}}} \right) }}{{{{\left[ {2BC\left( {\frac{\alpha }{C} - \rho - \frac{{\delta {\rho _1}}}{B}} \right) } \right] }^2}}} \\&- \frac{{2\left[ {{{(2BC)}^2}\theta \lambda b{\beta _1}\left( {\frac{2}{{{v^2}}}} \right) \left( {\frac{\alpha }{C} - \rho - \frac{{\delta {\rho _1}}}{B}} \right) {N_1}^\prime - {N_1}{{\left( {2BC\theta \lambda b{\beta _1}\left( {\frac{2}{{{v^2}}}} \right) } \right) }^2}} \right] }}{{{{\left[ {2BC\left( {\frac{\alpha }{C} - \rho - \frac{{\delta {\rho _1}}}{B}} \right) } \right] }^3}}} \\ \end{aligned} \end{aligned}$$

where

$$\begin{aligned} \begin{aligned} {N_1}&= C\left[ {2{{(\lambda b)}^2}{\beta _1}\left( {1 + \frac{{2\theta }}{v}} \right) + C{{(\lambda b)}^2}\left( {{\beta _2}{{\left( {1 + \frac{{2\theta }}{v}} \right) }^2} + \theta {\beta _1}\left( {\frac{2}{{{v^2}}}} \right) } \right) } \right] \\&\quad + \delta \left[ {2{{(\lambda b)}^2}{\beta _1}\left( {1 + \frac{{2\theta }}{v}} \right) + 2\lambda b{\rho _1} + C\left\{ \begin{aligned} {(\lambda b)^2}{\mu _2} + 2\lambda b{\beta _1}{\rho _1}\left( {1 + \frac{{2\theta }}{v}} \right) \\ \quad + {(\lambda b)^2}\left[ {{\beta _2}{{\left( {1 + \frac{{2\theta }}{v}} \right) }^2} + \theta {\beta _1}\left( {\frac{2}{{{v^2}}}} \right) } \right] \\ \end{aligned} \right\} } \right] , \\ {N_1}^\prime&= {(\lambda b)^2}BC\left[ {4\theta {\beta _2}\left( {\frac{{ - 1}}{{{v^2}}}} \right) \left( {1 + \frac{{2\theta }}{v}} \right) + 2\theta {\beta _1}\left( {\frac{{ - 2}}{{{v^3}}}} \right) } \right] \\&\quad + 4\theta {(\lambda b)^2}{\beta _1}B\left( {\frac{{ - 1}}{{{v^2}}}} \right) + 4\theta {\lambda b}{\beta _1}{\rho _1}\delta C\left( {\frac{{ - 1}}{{{v^2}}}} \right) , \\ {N_1}^{\prime \prime }&= {(\lambda b)^2}BC\left[ {4\theta {\beta _2}\left( {\frac{2}{{{v^3}}}} \right) \left( {1 + \frac{{2\theta }}{v}} \right) + 4\theta {\beta _1}\left( {\frac{3}{{{v^4}}}} \right) + 8{\theta ^2}{\beta _2}\left( {\frac{- 1}{v^2}}\right) ^2} \right] \\&\quad + 4\theta {(\lambda b)^2}{\beta _1}B\left( {\frac{2}{{{v^3}}}} \right) + 4\theta {(\lambda b)^2}{\beta _1}{\rho _1}\delta C\left( {\frac{2}{{{v^3}}}} \right) . \\ \end{aligned} \end{aligned}$$

Clearly all the partial derivatives of of TC exists and are continuous over the domain, thus the Hessian matrix of TC will be

$$\begin{aligned} M_1= \left[ \begin{array}{ll} \frac{{\partial ^2}{TC}}{{\partial \alpha }^2}&{} \frac{{\partial ^2}{TC}}{\partial \alpha \partial v} \\ \frac{{\partial ^2}{TC}}{\partial v \partial \alpha }&{} \frac{{\partial ^2}{TC}}{{\partial v}^2} \\ \end{array} \right] . \end{aligned}$$

Let the determinant of first and the second principle minor be denoted as

$$\begin{aligned} \begin{aligned} F_1&=\Big |\frac{{\partial ^2}{TC}}{{\partial \alpha }^2}\Big |,~~~~~ G_1=\Big |\left( \frac{{\partial ^2}{TC}}{{\partial \alpha }^2}\right) \left( \frac{{\partial ^2}{TC}}{{\partial v}^2}\right) -\left( \frac{{\partial ^2}{TC}}{\partial \alpha \partial v}\right) ^2\Big |\\ \end{aligned} \end{aligned}$$

respectively. As clearly from the 3D plots (Fig. 7a, b) we can see that \(F_1>0\) and \(G_1\geqq 0\) so \(M_1\) is a positive semidefinite matrix. Hence TC is a convex function. Thus (CP) is a convex programming problem over S. Therefore the optimal cost \(TC(\alpha ,v)=\$ {\mathbf {17}}\) obtained for the system is a global minima.

Fig. 7
figure 7

Surface plot of \(F_1\) and \(G_1\) in (a) and (b) respectively

Full size image

6.2 Bi-objective optimization

Most of the optimization studies of queueing system focus on a single-objective cost optimization problems. However, for real-world problems to make the system robust more than one objective function should be considered which will be conflicting in nature. Thus multiobjective progamming problem is of great importance. In this sub-segment, we formulate a bi-objective problem of optimization where the expected cost function (TC) and the expected awaiting time in queue \((W_1)\) are to be minimized simultaneously over a domain S. Consider the following bi-objective problem of optimization

$$\begin{aligned} (BP) ~~~~~~~~~~~~~~~~~~~~~(TC^*, {W_1}^*)= \min (TC, {W_1}) \end{aligned}$$
(19)

where \(S:\{(\alpha , v): \alpha , v\in (0.2,10)\}\subset R^2\).

In the problem (BP) there are two objective functions and thus we try to find the minimum cost to be invested with minimum awaiting time of the clients over the feasible region which is convex.

To attain the efficient solution of the problem (BP), the Multi-Objective Genetic Algorithm (MOGA) has been executed. It is to be noted from Pareto frontier that the Pareto optimal solution cannot improve both the objective functions simultaneously. Hence the concept of an efficient solution is studied. Figure 8 represents non-controlling solutions obtained using MOGA where it can be seen that for the Pareto rates \(\varvec{(\alpha ^*, v^*)=(0.9066, 9.9949)}\), the minimum cost and minimum client’s awaiting time are \(\mathbf {TC^*=\$17.4}\) and \(\mathbf {{W_1}^*=35.9}\) mins, respectively. Clearly, the minimum cost obtained for the bi-objective optimization problem (BP) has slightly increased from the minimum cost obtained for (CP) but as the minimum awaiting time of the clients is also obtained so, it will be very useful to develop a robust system.

Fig. 8
figure 8

Efficient solution using MOGA

Full size image

Table 2 highlights the changing pattern of the cost function and the halting time in the queue, simultaneously for changing values of \(\lambda \) and \(\theta \). It can be seen that as we increase the parametric values, the cost of the system increases while the halting time decreases which is in accordance with the Little’s formula and relates to the realistic scenario.

Table 2 Using MOGA, the effect of \(\lambda \) and \(\theta \) is obtained on system cost and waiting time
Full size table

To ensure that the minimum value obtained for (BP) is global Pareto optimal value, we need to check if the minimization problem (BP) is convex or not for which, both the objective functions need to be shown convex. As discussed in sub-segment. 6.1, the proposed objective TC is convex. To show the other proposed objective \(W_1\) convex, rewriting Eq. (13) we get

$$\begin{aligned} W_1(\alpha , v)=\frac{C_h}{\lambda }(Q_f+Q_g+Q_b+Q_c+Q_e). \end{aligned}$$
(20)

The second order partial derivatives of Eq. (20) will be

$$\begin{aligned} \begin{aligned} \frac{{\partial ^2}{W_1}}{{\partial \alpha }^2}&=\frac{C_h}{\lambda }\left( \frac{{\partial ^2}{Q_f}}{{\partial \alpha }^2}+\frac{{\partial ^2}{Q_g}}{{\partial \alpha }^2}+\frac{{\partial ^2}{Q_b}}{{\partial \alpha }^2}+\frac{{\partial ^2}{Q_c}}{{\partial \alpha }^2}+\frac{{\partial ^2}{Q_e}}{{\partial \alpha }^2}\right) ;\\ \frac{{\partial ^2}{W_1}}{\partial \alpha \partial v}&= \frac{C_h}{\lambda } \left( \frac{{\partial ^2}{Q_f}}{\partial \alpha \partial v}+\frac{{\partial ^2}{Q_g}}{\partial \alpha \partial v}+\frac{{\partial ^2}{Q_b}}{\partial \alpha \partial v}+\frac{{\partial ^2}{Q_c}}{\partial \alpha \partial v}+\frac{{\partial ^2}{Q_e}}{\partial \alpha \partial v}\right) ;\\ \frac{{\partial ^2}{W_1}}{{\partial v}^2}&= \frac{C_h}{\lambda } \left( \frac{{\partial ^2}{Q_f}}{{\partial v}^2}+\frac{{\partial ^2}{Q_g}}{{\partial v}^2}+\frac{{\partial ^2}{Q_b}}{{\partial v}^2}+\frac{{\partial ^2}{Q_c}}{{\partial v}^2}+\frac{{\partial ^2}{Q_e}}{{\partial v}^2}\right) ;\\ \end{aligned} \end{aligned}$$
(21)

where \(Q_f,Q_g, Q_b, Q_c, Q_e\) are as defined in sub-segment. 6.1. Clearly the second order partial derivatives of \(W_1\) exists and are continuous over the domain, thus the Hessian matrix is

$$\begin{aligned} M_2= \left[ \begin{array}{ll} \frac{{\partial ^2}{W_1}}{{\partial \alpha }^2}&{} \frac{{\partial ^2}{W_1}}{\partial \alpha \partial v} \\ \frac{{\partial ^2}{W_1}}{\partial v \partial \alpha }&{} \frac{{\partial ^2}{W_1}}{{\partial v}^2} \\ \end{array} \right] . \end{aligned}$$

Consider the determinant of the first principal minor as

$$\begin{aligned} F_2=\Bigg |\frac{{\partial ^2}{W_1}}{{\partial \alpha }^2}\Bigg | \end{aligned}$$

and the determinant of second principal minor as

$$\begin{aligned} G_2=\Bigg |\left( \frac{{\partial ^2}{W_1}}{{\partial \alpha }^2}\right) \left( \frac{{\partial ^2}{W_1}}{{\partial v}^2}\right) -\left( \frac{{\partial ^2}{W_1}}{\partial \alpha \partial v}\right) ^2\Bigg |. \end{aligned}$$

Figure 9a, b shows that \(F_2 > 0\) and \(G_2\geqq 0\). So \(M_2\) is a positive semidefinite matrix thus \(W_1\) is also a convex function. Therefore (BP) is a convex programming problem. Thus it can be deduced that the Pareto optimal values \(\mathbf {TC^*=\$17.4}, \mathbf {{W_1}^*=35.9}\) are global optimal values of a bi-objective optimization problem (BP).

Fig. 9
figure 9

Surface plot of first and second minor of \(W_1\) in (a) and (b) respectively

Full size image

It may be noted that in the optimization problems (CP) and (BP), the interval (0.2, 10) ensures the convexity of the objective functions and thus the problems are convex programming problems. If the lower bound or upper bound is changed, the Hessian matrix of the objective functions remain no longer positive semidefinite hence the convexity of the objective functions is disturbed which is in fact an important sufficiency condition.

7 Conclusion

In tackling various real life computing situations easily, cost-effective systems are made for which the optimal values of critical parameters of the system might provide an insight to network system designers, system operators and software system engineers. Application of the developed model is not only limited to mentioned in the paper rather it can be seen in call centres, mobile and communication workplaces or any organization using commuting systems, where erratic queueing system with retrial can play a crucial role in conveying quality help rapidly under decreased expense. The model M/G/1 having erratic server with delayed repair and balking could be adopted by computing managers in getting minimal cost with minimized awaiting time. The presumption of resisting behaviour of clients make the model flexible to work in real life. In real economical world, our model is reliable and more helpful as the minimum cost value obtained is not only found using heuristic technique but also the concept of convexity is applied to show the existence of global minima. In future, the model can be studied further by considering different types of vacations for the server. We can extend the paper in future by applying priority concept with multi optional services or one can think of using different control policies i.e. N-policy, F-policy on this model.