A Model of Optimal Production Planning Based on the Hysteretic Demand Curve

Abstract

The article considers a hysteretic model of consumer behaviour in mono-product markets. Demand generation with regard to an individual consumer is modeled using a non-ideal relay with inverted thresholds. Therefore, the sales rate is defined as an analogue of the Preisach converter. The article considers the problem of the optimal production, storage, and distribution of goods, taking into account the hysteretic nature of the demand curve. The problem is reduced to a non-classical optimal control problem with hysteretic non-linearities. The latter is solved using Pontryagin’s maximum principle. The adopted economic model is based on the binary relationship of consumers to the product: the product is bought or the product is not bought. Transitions between these states are determined within the framework of our model only by the price of the goods; therefore, only the operator of a non-ideal relay can accurately describe such a dependence. The article presents the results of computational experiments illustrating the theoretical assumptions.

1. Introduction

The production, storage, and distribution of goods is one of the most important problems in applied economics. A large number of production models use the demand curve, which characterizes the intention and the ability of consumers to purchase goods produced by a particular enterprise. The demand curve is most commonly defined as a deterministic dependence of the sales volume on the “instantaneous” parameters of the state of the market. Such dependencies usually allow for a qualitative description of the market dynamics over short periods of time. However, they do not take into account earlier values of the parameters and, thus, do not account for a number of important factors, including consumer inertia and demand shifts occurring over time [1].

1.1. Purpose/Background

A number of studies [2,3,4,5,6] have suggested the completion of dynamic models of socio-economic processes with hysteretic elements. Thus, ref. [5] argues for the use of hysteretic elements in dynamic models. Here, hysteresis relates to the sum of individual impacts of all the participants of the market. The effect of hysteretic non-linearities on the level of unemployment is studied in [7].

A recent study [8] investigated the robustness of a large variety of macroeconomic models. The authors demonstrated that non-linear operators, namely, the play operator and the Prandtl–Ishlinskii operator, can be used to design accurate models of financial and economic processes.

It is known that hysteresis is commonly observed in various fields of physics, chemistry, economics, and biology. It takes place when the state of a system is determined by external conditions at a particular moment in time, as well as at other moments in the past. In this case, separate parts of a complex system usually exhibit strongly non-linear properties, which affect the behaviour of the whole system. We should note that dynamics are rather difficult to describe due to the complexity of the phase space and non-differentiability of the corresponding operators. Nevertheless, there are a number of constructive models of hysteretic converters which have been used effectively in the modeling of technical systems [9,10,11,12]. The most recent studies of the problem are presented in [13,14].

A large number of studies are devoted to the problem of the optimal control of dynamic systems that have hysteresis blocks in their composition. Thus, in [15,16], the optimal control problem is solved for a system with two state variables, one of which has evolution regulated by a controlled ordinary differential equation, and the other contains a hysteresis operator (backlash, Prandtl–Ishlinski operator and Preisach operator). Using the dynamic programming method, the corresponding Hamilton–Jacobi equation of the first order is derived, and it is proved that the objective function is the only bounded uniformly continuous solution of the Cauchy problem.

In a number of studies, hysteresis terms are present both in the equations of state and of control [17,18]. Within the framework of the discrete-time approach, dynamic programming equations are obtained and a method for their numerical solution is presented. In the case of continuous time, the conclusion is essentially based on the concept of the derivatives of multivalued functions. These results can be applied to systems with sensors and actuators (for example, piezoceramic materials).

Some specific problems of mathematical physics, in particular, the diffusion problem, are considered in [19,20]. A constructive method based on a two-parameter penalty function is proposed. The small parameter determines the deviation of the solution in finite time from the expected value, while the second parameter is used to approximate the main variational inequalities. The solution to the problem of control can be obtained by taking the limit in a doubly degenerate control system; the convergence in spatial coordinates is observed in a strong sense and is uniform in time. Diffusion in biological systems is discussed in [21]. In this article, the analyzed system contains three diffusion equations describing the evolution of three biological species: prey, predator and prey food (or vegetation). The equation for food density includes a hysteresis operator in the form of a generalized stop. The problem of minimizing the integral functional of costs in relation to the solutions of the above-mentioned system is examined. Some relaxation-type results for the minimization problem are obtained and the existence of an almost optimal solution is established. The problem of reaction diffusion control [22] is modeled by a system with two types of control, namely distributed control functions and controls that act on part of the boundary of the domain. The equation of state is given by the reaction-diffusion system with the additional property that the scalar stop operator is included in the reaction rate. A feature of the output of the conjugate system is the non-locality in time of the Hadamard derivative of the state control operator.

Another interesting way to use hysteresis non-linearities in control problems is proposed in [23]. Here the hysteresis term is introduced as a variational inequality of evolution with a closed convex m-dimensional real arithmetic space. In this case, there are optimal solutions, as well as the necessary conditions for optimality of the first order. In particular, under certain assumptions about regularity, the behavior of solutions of the conjugate system is described in detail. It is emphasized that a significant difficulty in obtaining optimality conditions is caused by the non-differentiability of the non-linear operators considered.

In technical tasks, the control of elements of mechanical systems, in particular, the mechanical transmission, is of great interest [24,25]. In some cases, the task is reduced to eliminating the hysteresis effect in a high-precision mechanical transmission. To describe the hysteresis behavior of a harmonic drive, a modified Bouc–Wen model is used. With the help of coordinate transformation and feedback linearization, a mathematical model of a sequential mechanical drive system is obtained. The reference trajectory is tracked by a controller based on a linear quadratic controller. The relative error in tracking steady-state fluctuations tends to zero. When using optimal control, the output signal of the harmonic drive can follow a more complex trajectory. In a series of investigations, refs. [26,27,28,29] demonstrate the possibility of controlling magnetostrictive and piezoelectric drives. The main task is solved in detail, achieving a given trajectory of movement, as well as reducing vibrations.

In a recent paper [30], the problem of identifying the parameters of the Preisach converter was solved using the improved particle swarm optimization method, an important property of which is that it is not necessary to know the exact gradient of the optimized function. Using the example of calculating the characteristics of piezoelectric actuators (piezoelectric actuators), fast convergence, small calculation time and higher accuracy were demonstrated compared to the classical particle swarm optimization algorithm.

The Preisach model is one of the most prominent constructive models of hysteresis [31]. It is also one of the most popular and actively developing models, and employs hysteretic operators with stochastic parameters. Thus, [32] investigated the response of a non-linear system to stochastic external factors, while [33] modified the said model to the threshold numbers characterised by random, rather than deterministic variables. In [34], hysteretic operators with stochastic parameters are discussed.

Other hysteresis models, such as the Bouc–Wen model, can also be applied in this kind of approach. For example, in [26], the problem of controlling the drive of an atomic force microscope driven by piezoelectric actuators was solved. It was shown that the presence of hysteresis friction degraded the characteristics of the microscope and led to a loss of accuracy. To eliminate estimation errors, non-simulated vibration and disturbances, a sliding mode control with perturbation estimation was used. This method was used to improve the performance and reliability of the system—the designed controller can provide superior tracking performance for piezoelectric actuators with low input frequencies. Similar ideas have been considered in recent studies [35,36,37,38,39].

1.2. Method

To date, several methods for solving optimal control problems have been used: the Bellman principle, dynamic programming, and the Pontryagin maximum principle [40,41,42,43,44]. The first two of these require smoothness of the right parts of the dynamic system, as well as the integrand in the target functional. The problem considered in this paper contains strong non-linearities (i.e., non-linearities that do not allow the possibility of linearization). It is known that the operator corresponding to the Preisach converter does not have a weak Gâteaux derivative. Therefore, the Pontryagin maximum principle is used to solve the problem, where smoothness is not necessary [45].

It should be noted that a number of issues related to optimal control leading to non-classical variational problems have previously been considered by various authors. Further development of the theoretical foundations of control methods for complex systems has been aimed at the introduction of stochastic control and stochastic filtering of dynamic systems, and the construction of general methods for solving non-classical variational problems.

The adopted economic model is based on the binary relationship of consumers to the product: the product is bought or the product is not bought. Transitions between these states are determined within the framework of our model only by the price of the goods; therefore, only the operator of a non-ideal relay can accurately describe such a dependence.

1.3. Results and Conclusions

In this article, we consider the problem of the optimal production, storage, and distribution of goods, taking into account the hysteretic nature of the demand curve. Namely, the problem is considered with regard to the demand curve determined by means of the Preisach converter with an inverted individual relay. With certain limitations imposed on the parameters, the problem has a single solution.

To use the proposed model for the description of real economic processes, a method for identifying the model parameters is required. This requires statistical studies of the attitudes of individual consumers to certain goods; that is, each potential consumer needs to compare, by means of a survey, the threshold values of the price at which the goods either cease to be purchased or begin to be purchased. However, in this paper, we solve the problem associated with modeling hysteresis effects in economic relations, and also consider a relatively simple problem of optimizing the trade and production strategy of producers under the conditions of the hysteresis function of demand.

2. Non-Ideal Relay

To describe hysteresis non-linearity, special functional operators are introduced. This approach to describing the hysteretic behaviour of systems was suggested and actively developed by M. A. Krasnosel’skii and A. V. Pokrovskii [46]. The models of hysteretic operators are defined on the space of continuous functions, and the dynamics of these converters are described by the relations “input state” and “output state”.

Let $R[\alpha ,\beta ,x_0,t_0]$ denote a hysteretic converter corresponding to the non-ideal relay with thresholds $\alpha$ and $\beta$, where $x_0\in \{0,1\}$ is the initial state of the converter and $t_0$ is the initial moment. The state space of the non-ideal relay is a set of two elements $\{0,1\}$. The input of the system is a function $u\left(t\right)$, which is continuous, when $t⩾t_0$. The output of the system is a step function $x\left(t\right)$ determined by the operator relation

$$x\left(t\right)=R[\alpha ,\beta ,x_0,t_0]u\left(t\right).$$

(1)

The initial state $x_0$ of the converter has to comply with the following condition:

$$x_0=\left\{\begin{array}{c}0,\mathrm{if}u\left(0\right)⩽\alpha ,\hfill \\ 1,\mathrm{if}u\left(0\right)⩾\beta .\hfill \end{array}\right.$$

(2)

When $\alpha ⩽u\left(0\right)⩽\beta$ is satisfied, $x_0$ can take any value belonging to the set $\{0,1\}$. The output values $x\left(t\right)$, when the input $u\left(t\right)$ is continuous for $t\in (t_0,\infty )$ with any $t=\tau$, are determined according to the following rule:

$$R[\alpha ,\beta ,x_0,t_0]u\left(\tau \right)=\left\{\begin{array}{c}{x}_0,\mathrm{if}\forall t\in [t_0,\tau ]:[\alpha <u\left(t\right)<\beta ],\hfill \\ \phantom{-1}1,\mathrm{if}\exists t^{\prime }\in [t_0,\tau ):[u\left(t^{\prime }\right)⩾\beta ]\wedge \{\forall t\in [t^{\prime },\tau ]:[u\left(t\right)>\alpha ]\},\hfill \\ \phantom{-1}0,\mathrm{if}\exists t^{\prime }\in [t_0,\tau ):[u\left(t^{\prime }\right)⩽\alpha ]\wedge \{\forall t\in [t^{\prime },\tau ]:[u\left(t\right)<\beta ]\}.\hfill \end{array}\right.$$

(3)

Thus, the output takes the value of $[t_1,t_2]$, if either $x\left(t_1\right)=0$ and $u\left(t\right)⩽\alpha$ with $t\in [t_1,t_2]$, or $x\left(t_1\right)=1$ and $u\left(t\right)⩾\beta$ with $t\in [t_1,t_2]$. We can say that the relay is on, when the output is 1. Otherwise the relay is considered to be off. A detailed description of the non-ideal relay converter and its properties is given in [46].

The Preisach converter is a continuous analogue of non-ideal relays connected in parallel. The state space of the converter consists of pairs $\{u,z(\alpha ,\beta \left)\right\}$, where u is an arbitrary value and $z(\alpha ,\beta )$ is the characteristic function of the subset of the half-plane $\alpha <\beta$. For most applied problems, the support of the $z(\alpha ,\beta )$ function belongs to a finite set, although this is not a crucial limitation.

The output of this converter is determined by the relation:

$$R[z(\alpha ,\beta ),t_0,x_0]u\left(t\right)=\underset{\alpha <\beta }{\int \int }R[\alpha ,\beta ,x_0,t_0]u\left(t\right)\mathrm{d}\alpha \mathrm{d}\beta .$$

(4)

The state space of the converter is demonstrated in Figure 1.

3. Modeling the Demand Curve

As we have mentioned in the Introduction, there are a large number of models that describe—in more or less detail—the attitude of consumers to a particular good (a group of goods). These models do not always take into account the history of such relations. In this regard, we should note that [1,2,7] demonstrated the need to consider the ambiguity of economic processes. In this article, we suggest a consumer demand model which reflects the dependence of its present state on its history.

Let the sales rate $P\left(t\right)$ (i.e., the number of goods sold per unit of time) at any particular moment t depend on the price $c\left(t\right)$ of the goods and its previous values. We note that consumers’ attitude towards the goods may change over time. Let the attitude of a particular consumer be a function $\mathcal{R}\left[c\right(t\left)\right]$, which takes the values from the set $\{0,1\}$ according to the following rule:

$$\mathcal{R}\left[c\left(t\right)\right]=\left\{\begin{array}{c}1,\mathrm{if}c\left(t\right)⩽\alpha \left(t\right),\hfill \\ 0,\mathrm{if}c\left(t\right)⩾\beta \left(t\right),\hfill \\ 0\mathrm{or}1,\mathrm{if}\alpha \left(t\right)<c\left(t\right)<\beta \left(t\right),\hfill \end{array}\right.$$

(5)

According to (5), the function $\mathcal{R}\left[c\right(t\left)\right]$ equals 1, if the consumer purchases the goods at the moment t. Otherwise the function equals 0. Let us consider $\mathcal{R}\left[c\right(t\left)\right]$ to be the output of a converter $R[\alpha \left(t\right),\beta \left(t\right),R_0]$, analogous to the non-ideal relay with inverted thresholds $\alpha$, $\beta$, which receives an input signal $c\left(t\right)$ with $t\in [0,T]$, where $R_0$ is the initial state of the converter. The attitude towards the goods changing over time can be modeled as dependencies $\alpha =\alpha \left(t\right)$, $\beta =\beta \left(t\right)$.

The relations between the input and the output of the converter $R[\alpha ,\beta ,R_0]$ are schematically shown in Figure 2.

Let ${\gamma }_i$ denote the purchase rate of the i-th consumer, $i=1,\cdots ,n$. Then, for the system of n consumers, the sales function will be

$$\mathcal{P}\left(c\left(t\right)\right)=\sum _{i=1}^n{\gamma }_{i}R[{\alpha }_i,{\beta }_i,R_{0,\left(i\right)}]c\left(t\right).$$

(6)

Of practical importance is the case with a large number of individual consumers $n\to \infty$, when the impact of each participant of the market is ${\gamma }_i\to 0$. In this case, the sales function should be equal to the continuous limit of the sum (6):

$$\mathcal{P}\left(c\left(t\right)\right)=\underset{\alpha <\beta }{\int }\omega (\alpha ,\beta ,t)\mathrm{d}\mu ,$$

(7)

with designations

$$\omega (\alpha ,\beta ,t)=\mathsf{\Gamma }\left[{\omega }_0(\alpha ,\beta )\right]c\left(t\right)=R[\alpha \left(\gamma \right),\beta \left(\gamma \right),R_0\left(\gamma \right)]c\left(t\right),$$

(8)

and $\gamma \in \left\{\right(\alpha ,\beta ):\alpha <\beta \}$. The dependence of the measure $\mu =\mu \left(t\right)$ on the time makes it possible to update the dynamics of the attitude of groups of consumers towards the goods.

The converter (7) is analogous to the Preisach converter [46,47] with inverted 1 s and 0 s.

4. A Production Model

Below, we suggest that the support of the measure $\mu$ is contained in the triangle of the half-plane $\beta >\alpha$, described by the inequalities $\alpha >0$, $\beta <a$.

In the simplest case, the problem of profit maximisation at the current time is set as follows: it is necessary to determine the price function $c^{*}\left(t\right)$ on the time interval $[0,T]$ so that at the moment $t=t_1$, the profit

$$c\left(t\right)\mathcal{P}\left(t\right){|}_{t_1}\to max.$$

(9)

Let us assume that $\forall t\in [0,T]:0<c\left(t\right)⩽a$, i.e., the price of the goods is not negative and is limited from above by a.

This problem can be solved using the following function:

$$c^{*}\left(t\right)=\left\{\begin{array}{c}{\displaystyle \frac{a}{\tau }}\left(\tau -t\right),\mathrm{if}0⩽t<\tau ,\hfill \\ {\displaystyle \frac{a}{2(t_1-\tau )}}\left(t-\tau \right),\mathrm{if}\tau ⩽t⩽t_0.\hfill \end{array}\right.$$

(10)

The graph of the function $c^{*}\left(t\right)$ is shown in Figure 3a.

Indeed, at the moment $t=\tau$, the goods’ price is maximum $c\left(\tau \right)=a$, and all individual relays are off, which corresponds to the lack of demand. On the segment $[\tau ,t_1]$, the price monotonously decreases, and the demand $\mathcal{P}\left(c\left(t\right)\right)={\displaystyle \frac{a^2}{2}}-{\displaystyle \frac{c^2}{2}}$ grows (see Figure 3b). As we can see, when $t=t_1$, the profit reaches its maximum, since the derivative

$${\displaystyle \frac{\mathrm{d}\left(cP\right(c\left)\right)}{\mathrm{d}c}}={\displaystyle \frac{\mathrm{d}}{\mathrm{d}c}}\left(c\left({\displaystyle \frac{a^2}{2}}-{\displaystyle \frac{c^2}{2}}\right)\right)={\displaystyle \frac{1}{2}}\left(a^2-4c^2\right)$$

(11)

equals zero, when $c={\displaystyle \frac{a}{2}}$, and the second derivative ${\displaystyle \frac{{\mathrm{d}}^2\left(c\mathcal{P}\left(c\right)\right)}{\mathrm{d}{c}^2}}=-4c<0$, since the price is not negative.

Thus, for a rather simple production model, the profit maximisation is obtained through raising the price to its maximum level and then monotonously reducing it to $a/2$. However, a detailed analysis of a more realistic production model requires the use of optimal control methods [45,48].

Let us consider a problem of the optimal planning of the production process based on a hysteretic demand curve. For this, we should introduce the following designations (see Figure 4):

$Z\left(t\right)$—the number of goods at the manufacturer’s warehouse,

$V\left(t\right)$—the number of goods purchased by the consumer,

$U\left(t\right)$—the production rate, i.e., the number of goods supplied to the warehouse per unit of time,

$P\left(t\right)$—the sales rate, i.e., the number of goods sold per unit of time,

$k_1$—consumption coefficient,

$k_2$—coefficient characterising the cost of storage per goods item,

$c\left(t\right)$—price per goods item.

Let us assume that the production cost is 1. The system of integro-differential equations used to model the problem of production, storage, and distribution of goods is presented as follows:

$$\begin{array}{cc}\hfill & \dot{Z}=U-P,\hfill \end{array}$$

(12)

$$\begin{array}{cc}\hfill & \dot{V}=P-k_{1}V,\hfill \end{array}$$

(13)

$$\begin{array}{cc}\hfill & P\left(t\right)=Z\left(t\right)\underset{\alpha <\beta }{\int }\omega (\alpha ,\beta ,t)d\alpha d\beta ,\hfill \end{array}$$

(14)

$$\begin{array}{cc}\hfill & \omega (\alpha ,\beta ,t)=\mathsf{\Gamma }\left[{\omega }_0(\alpha ,\beta )\right]c\left(t\right),\hfill \end{array}$$

(15)

$$\begin{array}{cc}\hfill & Z\left(0\right)=V\left(0\right)=0.\hfill \end{array}$$

(16)

To solve the production optimisation problem, we need to introduce the following additional assumptions:

(1) $0⩽U\left(t\right)⩽U_0$ — the presence of the threshold $U_0$ of the maximum production rate;

(2) $0<c\left(t\right)⩽a$ — the price of the goods is not negative and is limited by a;

(3) the production cost equals 1 (by redetermining the unit of measurement of the number of goods, we can set the cost arbitrarily, provided that it is not negative).

Let us now consider the functional equal to the total income over a finite time interval T, taken with the opposite sign:

$$\mathcal{I}\left(T\right)={\int }_0^T[-c\left(t\right)P\left(t\right)+U\left(t\right)+k_{2}Z\left(t\right)]dt.$$

(17)

Now we need to determine the conditions which allow for the maximum functional $\mathcal{I}\left(T\right)$.

According to Pontryagin’s maximum principle [45,48], to reach the minimum of the functional (17), we need the maximum of the function $H(p_1,p_2,Z,V,P,U,c)$, called the Hamiltonian of system (16):

$$\begin{array}{c}H(p_1,p_2,Z,V,P,U,c)=p_1(U-P)+p_2(P-k_{1}V)-(-cP+U+k_{2}Z)=\hfill \\ \hfill =(p_1-1)U+(p_2-p_1+c)P-k_1{p}_{2}V-k_{2}Z.\end{array}$$

(18)

Since function H is linear along the argument U, $maxH(p_1,p_2,Z,V,P,U,c)=$

$H(p_1,p_2,Z,V,P,U^{*},c^{*})$ is reached when

$$U^{*}\left(t\right)=\left\{\begin{array}{c}0,\mathrm{if}{p}_1\left(t\right)<1,\hfill \\ U_0,\mathrm{if}{p}_1\left(t\right)>1,\hfill \\ \forall U,\mathrm{if}{p}_1\left(t\right)=1.\hfill \end{array}\right.$$

(19)

Let us describe the conditions, when summand $A\left(c\right)\equiv (p_2-p_1+c)P\left(c\right)$ makes an extreme contribution to the Hamiltonian (18). (1) Let $p_2-p_1+c>0$ be true. Then, $A\to max$ when $P=(a^2-c{\left(t\right)}^2)/2$ (see Figure 5a).

$${\displaystyle \frac{\mathrm{d}A\left(c\right)}{\mathrm{d}c}}={\displaystyle \frac{3}{2}}{c}^2-c(p_1-p_2)-{\displaystyle \frac{1}{2}}{a}^2.$$

(20)

The derivative equals zero when

$$c_{\pm }={\displaystyle \frac{1}{3}}[p_1-p_2\pm \sqrt{{(p_2-p_1)}^2+3a^2}].$$

(21)

According to the model, the price belongs to segment $[0,a]$. Therefore, we need only the positive root $c_{+}={\displaystyle \frac{1}{3}}[p_1-p_2+\sqrt{{(p_2-p_1)}^2+3a^2}]$. With all the possible values of $p_1,p_2,a$, this root belongs to segment $[0,a]$ and corresponds to the maximum of the function $A\left(c_{+}\right)$, since ${\displaystyle \frac{{\mathrm{d}}^{2}A\left(c\right)}{\mathrm{d}{c}^2}}{|}_{c=c_{+}}<0$.

(2) Let $p_2-p_1+c<0$ be true. Then, as shown in Figure 5b, $A\to max$ when $P={(a-c)}^2/2$.

$$A\left(t\right)=(p_2-p_1+c){\displaystyle \frac{{(a-c\left(t\right))}^2}{2}},$$

(22)

$${\displaystyle \frac{\mathrm{d}A\left(c\right)}{\mathrm{d}c}}={\displaystyle \frac{1}{2}}{(a-c)}^2-(a-c)(p_2-p_1+c).$$

(23)

The derivative equals zero when the price value equals either $c=a$, or

$$c={\displaystyle \frac{1}{3}}[a-2(p_2-p_1)].$$

(24)

Thus, the maximum of the Hamiltonian is reached at a certain $c^{*}\in [0,a]$, which satisfies the condition:

$$c^{*}=\left\{\begin{array}{c}{\displaystyle \frac{1}{3}}[p_1-p_2+\sqrt{{(p_2-p_1)}^2+3a^2}],\mathrm{if}{p}_2-p_1+c^{*}>0;\hfill \\ {\displaystyle \frac{1}{3}}[a-2(p_2-p_1)],\mathrm{if}{p}_2-p_1+c^{*}<0.\hfill \end{array}\right.$$

(25)

Taking into account the determined value of $c^{*}$, we obtain:

$$\underset{\{p_1,p_2,Z,V,P,U,c\}}{max}H=H^{*}(p_1,p_2,Z,V){|}_{c\left(t\right)=c^{*}}=(p_1-1)U^{*}+(p_2-p_1+c^{*})P-k_1{p}_{2}V-k_{2}Z.$$

(26)

According to Pontryagin’s maximum principle, the unknown functions $p_1\left(t\right)$ and $p_2\left(t\right)$ satisfy differential equations

$$\begin{array}{cc}\hfill & {\dot{p}}_1\left(t\right)=-{\displaystyle \frac{\partial H^{*}}{\partial Z}}=k_2-(p_2-p_1+c^{*}){\int }_{\alpha <\beta }\omega (\alpha ,\beta ,t)\mathrm{d}\alpha \mathrm{d}\beta ,\hfill \end{array}$$

(27)

$$\begin{array}{cc}\hfill & {\dot{p}}_2\left(t\right)=-{\displaystyle \frac{\partial H^{*}}{\partial V}}=k_1{p}_2\left(t\right)\hfill \end{array}$$

(28)

when $t\in [0,T]$ with boundary conditions $p_1\left(T\right)=p_2\left(T\right)=0$.

As can be seen, the solution to function $p_2\left(t\right)$ is trivial: $p_2\left(t\right)\equiv 0$ for any $t\in [0,T]$. Therefore, the dynamics of the integro-differential system (16) are determined by the evolution of the function $p_1\left(t\right)$ over time. The function, in turn, depends significantly on the output values of the hysteretic converter.

Relative to $p_1$, the equation is presented as:

$${\dot{p}}_1\left(t\right)=k_2-(c^{*}-p_1){\int }_{\alpha <\beta }\omega (\alpha ,\beta ,t)\mathrm{d}\alpha \mathrm{d}\beta .$$

(29)

Let us assume that the storage cost is low compared to the production cost and the revenues: $k_2=0$. In this case, the dependence of the parameter $p_1$ on the time in the interval $[0,T]$ is determined by the differential equation

$$\left\{\begin{array}{c}{\dot{p}}_1\left(t\right)=(p_1\left(t\right)-c^{*}\left(t\right))\underset{\alpha <\beta }{\int }\omega (\alpha ,\beta ,t)\mathrm{d}\alpha \mathrm{d}\beta ,\hfill \\ p_1\left(T\right)=0,\hfill \end{array}\right.$$

(30)

where

$$c^{*}=\left\{\begin{array}{c}{\displaystyle \frac{1}{3}}\left(p_1+\sqrt{p_1^2+3a^2}\right),\mathrm{if}{p}_1<c^{*},\hfill \\ {\displaystyle \frac{1}{3}}\left(a+2p_1\right),\mathrm{if}{p}_1>c^{*}.\hfill \end{array}\right.$$

(31)

According to (19), when $t=0$, the initial condition must satisfy the inequality $p_1\left(0\right)⩾1$. Indeed, we would otherwise observe the equality $U^{*}\left(0\right)=0$, which, according to the model, would mean that the production had not begun. Since the initial moment for our model is the start of the production process, the inequality $p_1\left(0\right)⩾1$ is true at the start of the production process.

The obtained inequality, in turn, imposes an important limitation on the relation (31): in the time interval $[0,T]$ only $p_1<c^{*}$ is implemented. Let us prove it.

Let us assume the opposite: $\exists t^{\prime }\in [0,T]:p_1\left(t^{\prime }\right)>c^{*}\left(t^{\prime }\right)$. Then, at this point $c^{*}\left(t^{\prime }\right)={\displaystyle \frac{1}{3}}\left(a+2p_1\right)$ or

$${\displaystyle \frac{1}{3}}(a-c\left(t^{\prime }\right))+{\displaystyle \frac{2}{3}}(p_1\left(t^{\prime }\right)-c\left(t^{\prime }\right))=0.$$

(32)

In (32), the second summand is positive, which means the first summand $(a-c\left(t^{\prime }\right))/3$$<0$. This contradicts the idea of price $c\left(t\right)$ in our model, since the price is limited over the whole observation range: $\forall t\in [0,T]:0<c\left(t\right)⩽a$. We can, thus, conclude that, for our model, only the first equation from (31) is implemented, namely:

$$c^{*}={\displaystyle \frac{1}{3}}\left(p_1+\sqrt{p_1^2+3a^2}\right).$$

(33)

Being inserted in (30) it results in:

$$\left\{\begin{array}{c}{\dot{p}}_1\left(t\right)={\displaystyle \frac{1}{3}}\left(2p_1\left(t\right)-\sqrt{p_1^2\left(t\right)+3a^2}\right){\int }_{\alpha <\beta }\omega (\alpha ,\beta ,t)\mathrm{d}\alpha \mathrm{d}\beta ,\hfill \\ p_1\left(T\right)=0,\hfill \\ 1<p_1\left(0\right)<a.\hfill \end{array}\right.$$

(34)

Therefore, we come to the following conclusion: in the production model (16) for any finite time moment T in the interval $[0,T]$, there is a single solution $c^{*}\left(t\right)$ to the system (34), which is strictly monotonous on $[0,T]$.

Indeed, we can easily see that $(\forall p_1<a):2p_1\left(t\right)-\sqrt{p_1^2\left(t\right)+3a^2}<0$. Therefore, extreme values of the functional $\mathcal{I}\left(T\right)$ are implemented by the function $p\left(t\right)$, for which ${\dot{p}}_1\left(t\right)<0$ and $p\left(T\right)=0$. It is obvious that this function monotonously decreases from a particular initial value $p_1\left(0\right)<a$ to zero. It follows from (31) that the value of $p_1$ on the right border corresponds to the price at the final moment: $c\left(T\right)=a/\sqrt{3}$.

This proves the existence and the strict monotony of the function $p_1\left(t\right)$ on $[0,T]$. Then, according to (33), the price is a strictly monotonous function of $p_1\left(t\right)$, which proves the existence of $c^{*}$ and its monotony in a strict sense. The uniqueness of $c^{*}$ stems from the uniqueness of the solution to the ordinary differential Equation (34) [49,50]. Therefore, for any final moment T system (34) has a single solution $c^{*}\left(t\right)$.

5. Results

We should point out that system (34) has a boundary condition corresponding to the end of the modeling time T. Therefore, the value of $p_1\left(0\right)$ is unknown. Unfortunately, the system (34) cannot be solved numerically by reversing the time axis due to the asymmetry of the Preisach converter with regard to the t parameter.

However, when $t=T$, the boundary condition can be satisfied using the following method: based on (34), we select the initial value $m\equiv p_1\left(0\right)\in [1,a]$ and perform the numerical solution. We can thus determine $p_1\left(T\right){|}_m$.

Let us consider the mapping $U^T:p^{\left(0\right)}\equiv p_1\left(0\right)\to p_1\left(T\right)$. Since $\exists p^{(-)}:(U^T{p}^{(-)}<0)$ and $\exists p^{(+)}:(U^T{p}^{(+)}>0)$, the continuity of the mapping $U^T$ results in $\exists p_1^{(*)}:(U^T{p}_1^{(*)}=0)$ [49,50].

Due to the monotony of $p_1\left(t\right)$ over the entire analysed segment, we can use a binary search to find the single root of the equation

$$p_1\left(T\right){|}_m=0.$$

(35)

The corresponding root m determines the initial value of the studied function at the moment $t=0$, matching the right border.

As we can see, the initial value of $p_1\left(0\right)$ depends on the parameters of the problem, namely on the state of the Preisach converter.

For the numerical experiments, we selected the following values of the dynamic system: $a=10$, $U_0=1$, $k_1=1.0$, $T\in \{10,50\}$.

Figure 6, Figure 7, Figure 8 and Figure 9 demonstrate the initial state of the hysteretic converter, corresponding to the demand at the moment $t=0$ and the function $p_1\left(t\right)$ and $c^{*}\left(t\right)$. The graphs show that, at the final moment, $p_1\left(t\right)$ turns zero and $c^{*}\left(T\right)$ takes the value $a/\sqrt{3}$. The graphs show the following functions: $U\left(t\right)$ (production rate), $Z\left(t\right)$ (the number of goods at the warehouse), $V\left(t\right)$ (the number of goods purchased by consumers), $P\left(t\right)$ (the sales rate).

The number of elementary relays in the analogue of the Preisach converter in our numerical calculations was $N=210$. The monotonous decrease in the price with the growth of t results in the monotonous increase in $\mathcal{P}\left(c\right(t\left)\right)$. Due to the finiteness of N, $\mathcal{P}\left(c\right(t\left)\right)$ changes rapidly at some points in time, which results in rapid changes in $P\left(t\right)=Z\left(t\right)\mathcal{P}\left(c\right(t\left)\right)$. The behaviour of $\mathcal{P}\left(c\right(t\left)\right)$ described above occurs after the “activation” of each of the subsequent lines of individual relays of the converter (see Figure 5b).

Figure 6, Figure 7, Figure 8 and Figure 9 demonstrate that, at the time of production ($\{t:U(t)>0\}$), the number of goods at the warehouse and the number of good purchased by consumers grow. When the production process stops, the derivative of the function $Z\left(t\right)$ becomes negative, while the growth of $V\left(t\right)$ continues. However, $V\left(t\right)$ soon reaches its maximum and begins to decrease monotonously due to the negative impact of the goods consumption $k_{1}V\left(t\right)$ in (12).

Let us also consider the value $Q\left(T\right)=-{\displaystyle \frac{\mathcal{I}\left(T\right)}{T}}$, showing the profit of the manufacturer per unit of time. The dependence of this value on the sales duration is presented in Figure 10. The monotonous dependence demonstrates that the relative profit grows with an increase in the duration of the production process.

6. Conclusions

The article proposes a model of demand generation based on the hysteretic approach. The attitude of consumers towards the goods was modeled using an operator analogous to the operator of a non-ideal relay with inverted thresholds. This model makes it possible to take into account the history of consumers’ attitudes at a finite time interval. The analogue of the said model for a finite number of consumers (when their number approaches infinity) is the Preisach converter, whose state space consists of characteristic functions of the half-plane $\{\alpha <\beta \}$ of a non-classical type and determines the clockwise direction of the hysteresis loop. In this article, we considered the problem of the optimal production, storage, and distribution of goods on mono-product markets, taking into account the hysteretic nature of the demand curve. The problem was reduced to a non-classical optimal control problem with hysteretic non-linearities. The latter was solved using the Pontryagin’s maximum principle. The study demonstrated that, with certain limitations imposed on the parameters, the problem has a single solution. The article presents model examples illustrating the search for an optimal solution to the studied problem.

Finally, we note that the results obtained can be used in solving optimization problems related to competition in single-commodity markets. It is also important to study open economic systems under the conditions of hysteresis behavior of economic agents.