SOME VARIATIONAL PROBLEMS RELATED TO THE LAUNCH OF AN ARTIFICIAL EARTH SATELLITE

D. E. Okhotsimsky, T. M. Eneev

This article considers the problem of placing an artificial Earth satellite into orbit. It is assumed that the launch is carried out by means of a rocket accelerator consisting of one or several stages. We investigate what the law of variation with time of the direction of thrust of the reactive engines should be in order to ensure the placement of the satellite into a prescribed orbit with the minimum expenditure of fuel. The most advantageous regime of fuel consumption is also sought.

The solution of the indicated problems is carried out under a number of simplifying assumptions and makes it possible to form a definite picture of the characteristic features of the optimal conditions for orbital insertion and to indicate ways along which one can achieve the construction of accelerators with minimum initial weight.

In § 1 the problem is considered of the simultaneous choice both of the program for the direction of thrust and of the fuel-consumption regime. In § 2 the problem of choosing the optimal program for a multistage accelerator with different numbers of stages is investigated under the assumption that the fuel-consumption regime is prescribed. In § 3 a generalization of the problem of placing a satellite into orbit is carried out for the case of motion in a central gravitational field with allowance for the rotation of the Earth.

§ 1. CHOICE OF THE OPTIMAL FUEL-CONSUMPTION REGIME AND THE OPTIMAL PROGRAM FOR THE DIRECTION OF THRUST

We shall solve the indicated problem under the assumption that aerodynamic forces are absent and that the field of terrestrial gravity is plane-parallel. The basis for the first hypothesis is the fact that, when placing a satellite into orbit, a significant part of the launch trajectory will lie in the upper layers of the atmosphere, where aerodynamic forces are small. Replacement of the Earth’s central field by a plane-parallel one is possible if it is assumed that the extent of the launch trajectory is small in comparison with the radius of the Earth.

Both of the indicated assumptions are in reality satisfied only approximately. However, the solution of the problem in the indicated approximate formulation is of undoubted interest, since it makes it very simple to obtain the complete solution of the problem, to analyze the result obtained, and to understand the basic regularities of the phenomenon.

The equations of motion of a satellite with thrusters, in projection onto rectangular coordinate axes, may be written in the form (Fig. 1)

\[ \begin{aligned} \frac{du}{dt} &= p\cos\varphi,\\ \frac{dw}{dt} &= p\sin\varphi - g,\\ \frac{dy}{dt} &= w,\\ \frac{dx}{dt} &= u, \end{aligned} \left. \right\} \tag{1.1} \]

where \(u\) and \(w\) are the horizontal and vertical projections of the velocity, \(p\) is the magnitude of the acceleration due to the reactive force, and \(\varphi\) is the angle of inclination of the thrust force to the horizon. Since the thrust may be regarded as directed along the longitudinal axis of the thruster, the angle \(\varphi\) may be understood as the angle of inclination of the thruster axis to the horizon (the pitch angle). In equations (1.1), \(x\) and \(y\) are the horizontal and vertical coordinates.

The last equation of system (1.1) serves only to determine the coordinate \(x\). If no restrictions are imposed on the flight range, then in solving the variational problem this equation may be omitted, since the coordinate \(x\) does not enter into the other equations. Let us consider the velocity \(V\) that the satellite would acquire in the ideal case, when no forces other than the reactive force acted on it. We have:

\[ V=\int_0^t p\,dt; \]

hence

\[ p=\frac{dV}{dt}. \]

Fig. 1.

Substituting into system (1.1), we write the first three equations in the form

\[ \begin{aligned} \frac{du}{dt} &= \frac{dV}{dt}\cos\varphi,\\ \frac{dw}{dt} &= \frac{dV}{dt}\sin\varphi - g,\\ \frac{dy}{dt} &= w. \end{aligned} \left. \right\} \tag{1.2} \]

The quantity \(V\) represents the available velocity that can be imparted to the satellite in the process of acceleration and insertion into orbit.

The right-hand sides of the equations of motion (1.2) contain two undetermined functions of time, \(V(t)\) and \(\varphi(t)\). We shall pose the problem of choosing these functions in such a way that, at the end of the insertion segment, the greatest horizontal velocity is obtained at a prescribed altitude. By virtue of reciprocity, the solution obtained will also ensure the attainment of the greatest altitude for a prescribed velocity, as well as the attainment of prescribed values of altitude and velocity with minimum fuel consumption.

Let us formulate the boundary conditions. Suppose that at the beginning of the motion, for \(t=0\), we have some altitude \(y_0\) and some values of the hori-

VARIATIONAL PROBLEMS CONNECTED WITH THE LAUNCHING OF AN EARTH SATELLITE

of the horizontal and vertical projections of the velocity, \(u_0\) and \(w_0\). The velocity \(V\) at the beginning of the motion is naturally taken to be zero. At the end of the motion, at the instant \(t=T\), the altitude must be equal to \(y=Y\), and the velocity must be directed horizontally, i.e. \(w=0\). The value of \(V\) at the end of the motion will be assumed equal to some fixed value \(V_k\). Thus it is supposed that there is a certain reserve of ideal velocity which must be used in the most rational way. If one assumes that the magnitude of the exhaust velocity, or of the specific thrust, does not depend on the fuel consumption per second, then the magnitude of the total ideal velocity \(V_k\) will be determined by the ratio of the initial and final weights of the accelerators and will not depend on the regime of fuel consumption. This means that prescribing a definite value of \(V_k\) at the end of the motion corresponds to prescribing a definite ratio of the initial and final weights, and, for a given satellite weight, means prescribing a definite initial weight of the accelerators.

Thus, the boundary conditions have the form:

\[ \begin{aligned} &\text{For } t=0 \qquad u=u_0,\quad w=w_0,\quad y=y_0,\quad V=0.\\ &\text{For } t=T \qquad \qquad\quad w=0,\quad y=Y,\quad V=V_k. \end{aligned} \tag{1,3} \]

The time of motion \(T\) on the injection segment may either be regarded as given, or it may be chosen from the condition of obtaining the greatest velocity at the end of the motion. In doing so one may assume that the time of motion \(T\) on the injection segment does not coincide with the operating time of the engine \(T_*\), but satisfies only the condition

\[ T \geq T_*, \]

where \(T_*\) is a prescribed quantity.

We shall not restrict the function \(\varphi(t)\) by any conditions except the condition of continuity. The function \(V(t)\), as the integral of a nonnegative quantity, should be regarded as nondecreasing. An approximate graph of the dependence of \(V\) on time is shown in Fig. 2. The horizontal portions correspond to motion with the engine not operating. The vertical portions correspond to instantaneous consumption of part of the fuel. If no restrictions are imposed on the magnitude of the fuel consumption per second, then any line in the \((t,V)\)-plane joining the points \(O\) and \(B\) (Fig. 2), along which \(V\) does not decrease, will be admissible. If restrictions are imposed on the consumption per second, then the admissible line will satisfy a more stringent condition, namely: at each of its points, the element of the curve issuing from that point must be situated inside or on the boundaries of a certain angle, the upper side of which corresponds to motion with the greatest, and the lower side with the least, consumption per second. In the limiting case, when the greatest admissible consumption is infinitely large (instantaneous burning), and the least is equal to zero (motion with the engine switched off), the upper side of the angle will go vertically, and the lower horizontally.

Fig. 2.

The velocity at the end of the acceleration segment is obtained by integrating the first equation (1,2) from \(0\) to \(T\):

\[ u=\int_0^T \frac{dV}{dt}\cos\varphi\,dt+u_0. \tag{1,4} \]

In addition, we have two other equations, which are differential constraints. We rewrite them in the form

\[ \frac{dw}{dt}-\frac{dV}{dt}\sin\varphi+g=0, \tag{1,5} \]

\[ \frac{dy}{dt}-w=0. \tag{1,6} \]

We also have the boundary conditions (1,3). We represent the quantities \(u_0\) and \(w_0\) in the form

\[ u_0=v_0\cos\vartheta_0;\qquad w_0=v_0\sin\vartheta_0. \tag{1,7} \]

Taking \(V_0\) as given and varying the angle \(\vartheta_0\), we obtain from the solution of the problem the value of the optimal angle of inclination of the initial-velocity vector.

Let us form the auxiliary functional

\[ J=v_0\cos\vartheta_0+\int_0^T\left\{\frac{dV}{dt}\cos\varphi+\lambda_1\left(\frac{dw}{dt}-\frac{dV}{dt}\sin\varphi+g\right)+\right. \]

\[ \left. +\lambda_2\left(\frac{dy}{dt}-w\right)\right\}dt, \tag{1,8} \]

where \(\lambda_1\) and \(\lambda_2\) are certain as yet undetermined functions of time. We compute the variation of the obtained functional, also varying the initial angle \(\vartheta_0\) and the time of motion on the launching segment \(T\).

Carrying out the variation with respect to all the functions and variable parameters entering it, and using the boundary conditions, we obtain the following expression for the variation:

\[ \delta J=-v_0(\sin\vartheta_0+\lambda_1\cos\vartheta_0)\delta\vartheta_0+g\lambda_1^{0}\delta T+ \]

\[ +\int_0^T\left\{\left(-\sin\varphi-\lambda_1\cos\varphi\right)\frac{dV}{dt}\delta\varphi+ \left[-\frac{d}{dt}\left(\cos\varphi-\lambda_1\sin\varphi\right)\right]\delta V-\right. \]

\[ \left. -\frac{d\lambda_2}{dt}\delta y-\left(\frac{d\lambda_1}{dt}+\lambda_2\right)\delta w\right\}dt . \tag{1,9} \]

If the initial angle is not fixed, but is selected from the optimum condition, then we have a relation between the initial angle \(\vartheta_0\) and the value of the function \(\lambda_1(t)\) at the beginning of the motion,

\[ \tg\vartheta_0=-\lambda_1\big|_{t=0}, \tag{1,10} \]

obtained by equating to zero the first of the expressions preceding the integral. Equating the second expression to zero, we obtain the condition at the end of the acceleration segment

\[ \lambda_1\big|_{t=T}=0. \tag{1,11} \]

If the time of motion is fixed, then \(\delta T=0\), and the condition written for the function \(\lambda_1(t)\) drops out.

We now determine \(\lambda_1\) and \(\lambda_2\) in such a way that the multipliers of \(\delta y\) and \(\delta w\) under the integral become zero and thereby these variations drop out of the expression for the variation of the functional. If this is done, then the expression for the variation takes the form

\[ \delta J=\int_0^T\left\{\left(-\sin\varphi-\lambda_1\cos\varphi\right)\frac{dV}{dt}\delta\varphi+ \left[-\frac{d}{dt}\left(\cos\varphi-\lambda_1\sin\varphi\right)\right]\delta V\right\}dt, \tag{1,12} \]

where the function \(\lambda_1\) must be determined from the equation

\[ \frac{d\lambda_1}{dt}=-\lambda_2;\qquad \frac{d\lambda_2}{dt}=0 \tag{1,13} \]

with the boundary condition (1,11), which should be taken into account in the case when \(T\) is not fixed.

From the second equation (1,13) we have

\[ \lambda_2=\mathrm{const}. \tag{1,14} \]

Substituting this value into the first equation (1,13), we obtain

\[ \lambda_1=C_1+C_2t, \tag{1,15} \]

where \(C_1\) and \(C_2\) are certain constants. In the case when \(T\) is not fixed, we have the relation

\[ C_1+C_2T=0 \]

and the expression for \(\lambda_1\) is reduced to the form

\[ \lambda_1=-C_2(T-t). \tag{1,16} \]

Formulas (1,15) and (1,16) show that the quantity \(\lambda_1\) is a linear function of time. Note that these formulas contain the same number of arbitrary constants, since in formula (1,16), in place of \(C_1\), the unknown quantity \(T\) enters.

Equating to zero the expression multiplying the variation \(\delta\varphi\) under the integral in formula (1,12), we obtain

\[ \operatorname{tg}\varphi=-\lambda_1 \tag{1,17} \]

and, using formula (1,15), we find the expression for the optimal program in tangage in the form

\[ \operatorname{tg}\varphi=-(C_1+C_2t). \tag{1,18} \]

If \(T\) is not fixed, then according to (1,16) we have

\[ \operatorname{tg}\varphi=C_2(T-t). \tag{1,19} \]

Consequently, if the time is not fixed, then at the end of the motion the tangage angle must be equal to zero, i.e. the axis of the accelerators must be directed horizontally.

The formula for the optimal program can also be rewritten in the form

\[ \operatorname{tg}\varphi=\operatorname{tg}\varphi_0-C_2t, \tag{1,20} \]

where \(\varphi_0\) is the tangage angle at the initial instant of time.

Thus, we arrive at the conclusion that, in the optimal program, the tangent of the tangage angle must be a linear function of time. The two parameters entering into the general expression for the program must be chosen so as to satisfy the boundary conditions.

Comparison of formulas (1,10) and (1,17) shows that the optimal initial angle \(\vartheta_0\) must satisfy the condition

\[ \operatorname{tg}\vartheta_0=\operatorname{tg}\varphi_0. \tag{1,21} \]

This means that the direction of the longitudinal axis of the accelerator at the initial instant must coincide with the direction of the initial velocity. In the case where the angle \(\vartheta_0\) is prescribed in advance, this condition may fail to be fulfilled.

Formulas (1,18)—(1,20) give the law of variation of the tangage angle for any dependence \(V(t)\). We shall now choose the function \(V(t)\)

in the best way, assuming that for each such function the best law of variation of the pitch angle has been chosen each time.

Then the variation of the functional will take the form

\[ \delta J=\int_{0}^{T}\left\{-\left[\frac{d}{dt}\left(\sin\varphi-\lambda_{1}\cos\varphi\right)\right]\delta V\right\}\,dt . \tag{1.22} \]

When the differential conditions (1.5) and (1.6) are fulfilled, the functionals (1.4) and (1.8) coincide. If these conditions are identically satisfied under variation, then the variations of the functionals must also coincide. Therefore formula (1.22) makes it possible to calculate the variation of the velocity \(\delta U\) at the end of the acceleration segment when the fuel-consumption regime and the associated function \(V(t)\) are changed. It should be noted that the expression for the variation (1.22) will be valid only in the case when the entire change in \(V\) ends before the instant \(T\). Formula (1.22) can be transformed to the form

\[ \delta U=\int_{0}^{T}\Phi(t)\cdot\delta V\cdot dt, \tag{1.23} \]

where

\[ \Phi=\frac{C_{2}(\operatorname{tg}\varphi_{0}-C_{2}t)}{\sqrt{(\operatorname{tg}\varphi_{0}-C_{2}t)^{2}+1}} \tag{1.24} \]

or

\[ \Phi=C_{2}\sin\varphi. \tag{1.24a} \]

Fig. 3.

Formulas (1.23) and (1.24) show that the function \(\Phi(t)\), which stands under the integral as a multiplier of the variation \(\delta V\), depends only on time. On the interval \((0,T)\) the function may either have a constant sign throughout, or change sign; moreover, formula (1.24) shows that the change of sign of \(\Phi\) is associated with a change of sign of the numerator and therefore can occur only once.

On the basis of the obtained properties of the function \(\Phi(t)\), it is easy to obtain that the extremum of the functional is attained on functions \(V(t)\) with jump discontinuities. Let us consider separately the various possible cases.

a) The expression \(\Phi\) is positive throughout (Fig. 3, a). In this case the burning of the entire fuel reserve must occur instantaneously at the very beginning of the motion, since any other admissible line in the plane \((t,V)\) can otherwise be varied upward, as a result of which we obtain \(\delta U>0\), which means an increase of the velocity at the end of the ascent segment.

b) The expression \(\Phi\) is negative throughout (Fig. 3, b). Analogously to the preceding case, we obtain that in this case the best type of fuel consumption will be the instantaneous expenditure of the entire fuel reserve at the very end of the ascent segment, at \(t=T\).

c) The sign of \(\Phi\) changes from plus to minus (Fig. 3, c). In this case the extremum can be attained only on a broken line consisting of vertical, horizontal, and vertical segments. Any other admissible line can be varied in each of the regions of constancy of sign in such a way that the increment of the final velocity will be positive. In order to determine the position of the horizontal segment, let us write the expression for the variation of the final velocity under

vertical displacement of the segment. We obtain

\[ \delta U=\delta V\int_{0}^{T}\Phi(t)\,dt. \tag{1.25} \]

If the value of the integral is positive, then this case is analogous to case a); if it is negative, to case b). If the value of the integral is zero, then the extremum of the functional may be attained for some intermediate position of the horizontal segment. This means that, under certain conditions, it may prove advantageous to expend part of the fuel instantaneously at the very beginning of the motion, and the remaining part instantaneously at the very end of the motion.

d) The sign of \(\Phi\) changes from minus to plus (Fig. 3, d). In this case the expenditure of the entire fuel reserve should be made instantaneously at some moment of the motion.

All the possibilities indicated may, under certain conditions, be realized.

Let us consider the case when the time of motion is not fixed, but is chosen from the condition of maximum final velocity. The initial angle is likewise not fixed and is chosen optimally. In this case, according to (1.11), (1.17), and (1.20), we have

\[ \tg\varphi_k=\tg\varphi_0-C_2T=0, \tag{1.26} \]

whence

\[ C_2=\frac{\tg\varphi_0}{T}. \tag{1.27} \]

In this case the numerator in formula (1.24) for \(\Phi\) will have the form

\[ \tg^2\varphi_0\left(1-\frac{t}{T}\right) \tag{1.28} \]

and will be positive over the entire interval \((0,T)\). This means that case a) will be realized. We see that the greatest value of the final velocity \(U\), in the case when the initial angle and the time of motion are chosen from the optimum condition, is attained when the entire fuel reserve is expended instantaneously at the very beginning of the motion. The directions of the initial velocity \(v_0\) and of the additional velocity \(V\) must coincide. Here and in what follows we omit the subscript \(k\) on the quantity \(V\).

In what follows we shall assume that the initial velocity \(v_0\) is included in \(V\), and that the division of the total velocity reserve into two parts, if this is needed, is carried out in the optimal manner. It is then evidently possible to restrict ourselves to considering the case \(v_0=0\). We shall investigate this case both for free and for fixed time of motion on the launch segment. We shall determine for what relations of the parameters a solution is possible, and investigate the dependence of the character of the optimal motion on the magnitude of the prescribed launch time \(T\).

It is clear that for \(v_0=0\) only cases a and c are possible. To satisfy the boundary conditions we have the relations:

\[ V_1\sin\varphi_1+V_2\sin\varphi_2-gT=0, \tag{1.29} \]

\[ V_1\sin\varphi_1\cdot T-g\frac{T^2}{2}=Y, \tag{1.30} \]

where \(Y\) is the prescribed height, \(V_1\) is the velocity obtained at the initial moment of motion, and \(V_2\) is the velocity obtained at the final moment of motion. The sum

\[ V_1+V_2=V, \tag{1.31} \]

where \(V\) is the prescribed reserve of available velocity.

If \(T\) is not fixed, then \(V_2=0\). Equations (1.29) and (1.30) give the relation

\[ \frac{gT^2}{2}=Y, \tag{1.32} \]

which serves to determine the time \(T\):

\[ T=\sqrt{\frac{2Y}{g}}. \tag{1.33} \]

Substituting this value into either of equations (1.29) or (1.30), we obtain

\[ \sin\varphi_1=\frac{\sqrt{2gY}}{V}. \tag{1.34} \]

The problem is solvable if the right-hand side is less than unity. This means that, for a given available velocity \(V\), the altitude of injection \(Y\) must not be taken too large and, conversely, for a given altitude the reserve of available velocity must be sufficiently large; in any case, it must be greater than the velocity required to throw the body vertically upward to the height \(Y\). When this condition is fulfilled, the horizontal velocity at the end of the injection segment will be

\[ U=V\cdot \cos\varphi_1, \tag{1.35} \]

or

\[ U=V\sqrt{1-\frac{2gY}{V^2}}. \tag{1.36} \]

If \(T\) is fixed, then condition (1.32), generally speaking, is not satisfied, and no solution can be found for \(V_2=0\). The solution must be sought according to scheme c), with an intermediate position of the horizontal segment.

It is easy to see that the expression for the function \(\Phi\) can be represented in the form

\[ \Phi=\frac{d}{dt}\left(\frac{1}{\cos\varphi}\right). \tag{1.37} \]

Therefore the equality to zero of the integral in formula (1.25) gives

\[ \cos\varphi_1=\cos\varphi_2, \tag{1.38} \]

whence

\[ \varphi_2=\pm\varphi_1. \tag{1.39} \]

Let us consider both possibilities and determine under what relations between the parameters of the problem they are realized.

1. Let \(\varphi_2=\varphi_1\). From formula (1.29) we have

\[ \sin\varphi_1=\frac{gT}{V}, \tag{1.40} \]

and on the basis of (1.30) and (1.31) we obtain

\[ V_1=\frac{V}{2}\left(1+\frac{2Y}{gT^2}\right), \tag{1.41} \]

\[ V_2=\frac{V}{2}\left(1-\frac{2Y}{gT^2}\right). \tag{1.42} \]

Since \(V_2\) must be positive, we see that the given case occurs if the prescribed time of motion on the injection segment satisfies the condition

\[ T>\sqrt{\frac{2Y}{g}}, \tag{1.43} \]

i.e., greater than the optimal insertion time. For the final velocity we have the formula

\[ U=V\sqrt{1-\left(\frac{gT}{V}\right)^2}. \tag{1.44} \]

The condition that the expression \(\dfrac{gT}{V}\) must be less than unity is a requirement imposed on the magnitude of the available velocity reserve.

1. Let \(\varphi_2=-\varphi_1\). Then we have

\[ \sin \varphi_1=\frac{2Y}{VT}. \tag{1.45} \]

For the velocities we obtain

\[ V_1=\frac{V}{2}\left(1+\frac{gT^2}{2Y}\right), \tag{1.46} \]

\[ V_2=\frac{V}{2}\left(1-\frac{gT^2}{2Y}\right). \tag{1.47} \]

The solution exists if the prescribed time \(T\) is less than the optimal one, i.e., if

\[ T<\sqrt{\frac{2Y}{g}}. \tag{1.48} \]

The final velocity is computed by the formula

\[ U=V\sqrt{1-\left(\frac{2Y}{VT}\right)^2}. \tag{1.49} \]

Fig. 4.

The condition that the expression \(\dfrac{2Y}{VT}\) must be less than unity, analogously to the preceding case, may be regarded as a requirement that the available velocity reserve \(V\) be sufficiently large.

The formulas obtained show that the magnitude of the additional velocity \(V_2\), imparted at the end of the motion on the insertion segment, depends on how much the prescribed insertion time differs from the optimal one, and will be the larger, the greater this difference. If the prescribed time is equal to the optimal one, then \(V_2=0\), and both cases reduce to the preceding one. Let us note that in both cases

\[ V_2<V_1, \]

i.e., the second additional velocity is always smaller in magnitude than the first.

Conditions (1.43) and (1.48) mean that the velocity increment must occur, in the first case, on the descending branch, and in the second case on the ascending branch of the parabola obtained in the motion with initial velocity \(V_1\). Both cases are shown in Fig. 4. When the prescribed time is greater than the optimal one, the velocity increment \(V_2\) is directed parallel to the initial velocity \(V_1\). When, however, the prescribed insertion time is less than the optimal one, the velocity increment \(V_2\) is directed downward at an angle to the horizontal equal in magnitude to the initial angle \(\varphi_1\).

On the insertion segment the problems of ascent to a prescribed height and imparting the necessary velocity in the horizontal direction must be solved. The results obtained show that, in the most

in the advantageous case when the injection time has also been chosen optimally, the acquisition of velocity and the acquisition of altitude should be carried out by applying an instantaneous impulse at the very beginning of the motion. Injection into the horizontal direction of motion must be accomplished by the action of gravity. If the prescribed injection time differs from the optimum, then for injection into the horizontal direction it is necessary to use an impulse imparted at the end of the injection segment. If the prescribed time is excessively large, then the deflection of the velocity vector due to gravity is greater than necessary, and the vertical component of the second impulse must be directed upward. If the prescribed time is small, gravity does not have time to turn the velocity vector, and the deficient turn must be produced by means of the second impulse, which has a vertical component directed downward.

In obtaining the solution we proceeded from the assumption that an impulse may be imparted instantaneously and that the entire fuel supply, or part of it, may be expended instantaneously. Since in reality this is impossible, the best solution of the stated injection problem will be obtained by replacing instantaneous combustion by combustion of fuel at the maximum fuel consumption per second, and the interval between the application of the impulses by motion with the minimum fuel consumption per second or, better still, by motion with the engines switched off.

It was obtained above that, for motion in a plane-parallel gravitational field, the best injection is accomplished by applying one impulse at the beginning of the motion. For injection into an orbit in a central field of forces, such a solution proves unacceptable, since in this way it is impossible to obtain, for example, a circular orbit or an orbit that neither intersects nor is tangent to the surface of the Earth. Therefore, for applying the results obtained above to the solution of the problem of injection in a central field of forces, the problem of injection into orbit in a prescribed time is of fundamental importance. If the prescribed time is not too large, then the dimensions of the injection trajectory will be small in comparison with the radius of the Earth, and the hypotheses underlying the approximate treatment will not be strongly violated.

Let us compute an example which makes it possible to estimate the magnitude of the extent of the injection trajectory and the magnitude of the required supply of disposable velocity for various values of the prescribed injection time. In computing the example we shall restrict ourselves to the most important case, when the prescribed injection time does not exceed the optimum injection time in a plane-parallel field of forces.

Suppose it is necessary to inject a satellite to an altitude \(Y = 300\) km with horizontal velocity \(U = 7900\) m/sec.

By the formula

\[ T=\sqrt{\frac{2Y}{g}} \]

we obtain the optimum injection time in a plane-parallel field of forces

\[ T_{\mathrm{opt}} = 248\ \text{sec}. \]

In Fig. 5, for the region \(0.4\,T_{\mathrm{opt}} < T < T_{\mathrm{opt}}\), a graph is given of the dependence on \(T\) of the total required disposable velocity \(V\), the initial impulse \(v_1\), the initial angle \(\varphi_1\), and the horizontal range \(X\). From the figure it is evident that, as the injection time is shortened, the total required velocity increases, the magnitude of the initial impulse decreases, and that of the final impulse increases. The magnitude of the initial angle increases as the injection time decreases; the trajectories become increasingly more

steep, and the horizontal range of the insertion segment decreases. The obtained values of the horizontal range are considerably less than the radius of the Earth. Therefore there is sufficient reason to regard the example presented as an approximate calculation of the optimal insertion of a satellite into orbit.

It should be noted, in conclusion of the paragraph, that, as was shown above, the linear law of variation of the tangent of the pitch angle with time does not depend on the character of the variation of the function \(V(t)\). It will be applicable when using both single-stage and multistage accelerators, and also in the case of intervals between the end of operation of one stage and the beginning of operation of another. This result will be used in the next paragraph when considering insertion by means of an accelerator of the staged-rocket type.

Fig. 5.

It should also be noted that, in other formulations of the variational problem, the law for the tangent of the pitch angle may turn out to differ from a linear one. Thus, for example, if, in choosing the optimal time of motion, the horizontal range on the insertion segment is fixed, then the optimal law of variation in time of the tangent of the pitch angle turns out to be piecewise linear.

§ 2. MOTION WITH A GIVEN FUEL-CONSUMPTION REGIME

It was obtained above that the best control law for pitch is ensured by a linear dependence of the tangent of the pitch angle on time

\[ \operatorname{tg}\varphi = - (C_1 + C_2 t). \tag{2,1} \]

Integrating the equations of motion (1,1)—(1,4), we obtain expressions for the velocity projections and for the coordinates at some instant of time \(t_1\):

\[ u = \int_0^{t_1} p(t)\cdot \cos\varphi \cdot dt + u_0, \tag{2,2} \]

\[ w = \int_0^{t_1} p(t)\cdot \sin\varphi \cdot dt + w_0 - gt_1, \tag{2,3} \]

\[ x = \int_0^{t_1} (t_1 - t)\cdot p(t)\cdot \cos\varphi \cdot dt + u_0 t_1 + x_0, \tag{2,4} \]

\[ y = \int_0^{t_1} (t_1 - t)\cdot p(t)\cdot \sin\varphi \cdot dt + w_0 t_1 + y_0 - \frac{g t_1^2}{2}. \tag{2,5} \]

The quantities \(\cos\varphi\) and \(\sin\varphi\) may, on the basis of (2.1), be represented as functions of time

\[ \cos\varphi=\frac{1}{\sqrt{1+(C_1+C_2t)^2}}, \tag{2.6} \]

\[ \sin\varphi=-\frac{C_1+C_2t}{\sqrt{1+(C_1+C_2t)^2}}. \tag{2.7} \]

Since the function \(p(t)\) is regarded as given, the integrals in the right-hand sides of formulas (2.2)—(2.5) can, generally speaking, be computed.

The right-hand sides of formulas (2.2)—(2.5) contain two arbitrary constants associated with the control program. By choosing these constants, one can, generally speaking, ensure that, for a fixed time of motion along the acceleration segment, the specified altitude \(y_k\) and equality to zero of the vertical component of velocity will be attained at the end of the acceleration segment. Having thus chosen the constants, one can then determine, from formulas (2.2) and (2.4), the magnitude of the velocity \(u_k\) and the quantity \(x_k\).

For the actual computation of the integrals in the right-hand sides of formulas (2.2)—(2.5), it is necessary to know the dependence of the reactive acceleration on time, i.e., to know the function \(p(t)\). We note that the function \(p(t)\) is in general discontinuous, for example in the case of a composite rocket. Since composite rockets for an artificial satellite are of greatest interest, a more detailed investigation of the optimal motion of a composite rocket will be carried out below.

For the case of the motion of a composite rocket, using (2.2)—(2.5), we shall represent the kinematic parameters at the end of the injection trajectory in the form

\[ u_k=u_0+\sum_{i=1}^{n}u_k^{(i)}, \tag{2.8} \]

\[ w_k=w_0+\sum_{i=1}^{n}w_k^{(i)}, \tag{2.9} \]

\[ x_k=x_0+\sum_{i=1}^{n}x_k^{(i)}, \tag{2.10} \]

\[ y_k=y_0+\sum_{i=1}^{n}y_k^{(i)}, \tag{2.11} \]

where

\[ u_k^{(i)}=\int_{t_0^{(i)}}^{t_k^{(i)}} p_i(t)\cos\varphi\,dt, \tag{2.12} \]

\[ w_k^{(i)}=\int_{t_0^{(i)}}^{t_k^{(i)}} p_i(t)\sin\varphi\,dt -g\bigl(t_k^{(i)}-t_0^{(i)}\bigr), \tag{2.13} \]

\[ x_k^{(i)}=\int_{t_0^{(i)}}^{t_k^{(i)}}\bigl(t_k^{(i)}-t\bigr)p_i(t)\cos\varphi\,dt +\bigl(t_k^{(i)}-t_0^{(i)}\bigr)\cdot\sum_{q=0}^{i-1}u_k^{(q)}. \tag{2.14} \]

$$ y_k^{(i)}=\int_{t_0^{(i)}}^{t_k^{(i)}}\left(t_k^{(i)}-t\right)\cdot p_i(t)\cdot\sin\varphi\cdot dt+\left(t_k^{(i)}-t_0^{(i)}\right)\times $$

$$ \times\sum_{q=0}^{i-1} w_k^{(q)}-\frac{1}{2}g\left(t_k^{(i)}-t_0^{(i)}\right), \tag{2.15} $$

$$ u_k^0=u_0,\qquad w_k^0=w_0. $$

Here \(t_0^{(i)}\) and \(t_k^{(i)}\) are the times at the initial and final instants of motion, \(u_k^{(i)}\) and \(w_k^{(i)}\) are the increments of the horizontal and vertical components of velocity, \(x_k^{(i)}\) and \(y_k^{(i)}\) are the increments of the horizontal and vertical coordinates during the motion, and \(p_i(t)\) is the acceleration due to the reactive force; all quantities refer to the motion of the \(i\)-th stage of the rocket.

It is obvious that formulas (2.12)–(2.15) also cover the case when pauses occur between the operation of the engines of the separate stages of the composite rocket. In this case one may assume that the time interval corresponding to the pause belongs to a certain stage of the rocket, for which

$$ p_i(t)\equiv 0. $$

Let us carry out the calculation of the integrals in the right-hand sides of formulas (2.12)–(2.15) for one very important special case, in which the motion of each stage of the rocket occurs with constant thrust. In this case, for the \(i\)-th stage we have

$$ p_i(t)=\frac{gP_i}{G_i}, \tag{2.16} $$

where \(P_i\) is the engine thrust and \(G_i\) is the weight. The weight \(G_i\) can be represented in the form

$$ G_i=G_0^{(i)}\cdot\mu^{(i)} =G_0^{(i)}\left[1-\left(1-\mu_k^{(i)}\right)\cdot\frac{t-t_0^{(i)}}{T_i}\right], \tag{2.17} $$

where \(G_0^{(i)}\) is the initial weight, \(\mu_k^{(i)}\) is the final mass ratio, and \(T_i\) is the engine operating time. Introducing the notation

$$ \nu_0^{(i)}=\frac{G_0^{(i)}}{P_i}, \tag{2.18} $$

we obtain the expression for \(p_i(t)\)

$$ p_i(t)=\frac{g}{\nu_0^{(i)}\cdot\left[1-\left(1-\mu_k^{(i)}\right)\frac{t-t_0^{(i)}}{T_i}\right]}. \tag{2.19} $$

Determining the engine thrust by the formula

$$ P_i=-\frac{c_i}{g}\cdot\frac{dG_i}{dt}, \tag{2.20} $$

where \(c_i\) is the velocity of the outflow of gases from the engine nozzle, we easily obtain for the time \(T_i\) the formula

$$ T_i=\frac{c_i}{g}\nu_0^{(i)}\left(1-\mu_k^{(i)}\right). \tag{2.21} $$

We note that, for constant engine thrust, the mass flow per second is constant and the rocket mass is a linear function of time. For each stage one may therefore regard the tangent of the pitch angle as depending linearly not on time, but on the mass of the stage or on its relative mass. We then have

\[ \tg \varphi = A_i + B_i \mu^{(i)}, \tag{2.22} \]

where \(A_i\) and \(B_i\) are certain constants. At the beginning of the motion of the stage \(\mu^{(i)}=1\), and

\[ \tg \varphi_0^{(i)} = A_i + B_i . \tag{2.23} \]

At the instant when the operation of the stage engine ends, \(\mu^{(i)}=\mu_k^{(i)}\), and

\[ \tg \varphi_k^{(i)} = A_i + B_i \mu_k^{(i)} . \tag{2.24} \]

Here \(\varphi_0^{(i)}\) and \(\varphi_k^{(i)}\) are the values of the pitch angle at the beginning and at the end of the motion of the stage.

Substituting into formulas (2.8)—(2.15) the expression for \(p_i(t)\) from (2.19) and the expressions for \(\sin \varphi\) and \(\cos \varphi\) from (2.6) and (2.7), we obtain after integration

\[ u_k = u_0 + \sum_{i=1}^{n} c_i \cdot f_1^{(i)}, \tag{2.25} \]

\[ w_k = w_0 + \sum_{i=1}^{n} c_i \cdot \left[f_2^{(i)} - v_0^{(i)} \cdot \left(1-\mu_k^{(i)}\right)\right], \tag{2.26} \]

\[ x_k = x_0 + \sum_{i=1}^{n} \frac{c_i^2}{g}\, v_0^{(i)} \cdot \left[ f_3^{(i)} + \left(1-\mu_k^{(i)}\right) \cdot \sum_{q=0}^{i-1} \frac{c_q}{c_i}\, f_1^{(q)} \right], \tag{2.27} \]

\[ \begin{aligned} y_k = y_0 + \sum_{i=1}^{n} \frac{c_i^2}{g}\, v_0^{(i)} \Bigg\{ & f_4^{(i)} -\frac{1}{2} v_0^{(i)} \left(1-\mu_k^{(i)}\right) \\ &+\left(1-\mu_k^{(i)}\right) \cdot \sum_{q=0}^{i-1} \frac{c_q}{c_i} \cdot \left[ f_2^{(q)} - v_0^{(q)} \cdot \left(1-\mu_k^{(q)}\right) \right] \Bigg\}, \end{aligned} \tag{2.28} \]

where the dimensionless functions \(f_1^{(i)}, f_2^{(i)}, f_3^{(i)}\), and \(f_4^{(i)}\) are determined by the expressions

\[ f_1^{(0)}=\frac{u_0}{c_0}, \qquad f_2^{(0)}=\frac{w_0}{c_0}, \]

\[ \begin{aligned} f_1^{(i)} = \frac{1}{\sqrt{1+A_i^2}} \Bigg\{ &\operatorname{Arcsh} \left( A_i+\frac{1+A_i^2}{\mu_k^{(i)} B_i} \right) \\ &- \operatorname{Arcsh} \left( A_i+\frac{1+A_i^2}{B_i} \right) \Bigg\}, \end{aligned} \tag{2.29} \]

\[ f_2^{(i)} = A_i f_1^{(i)} + \operatorname{Arcsh}(A_i+B_i) - \operatorname{Arcsh}\left(A_i+B_i\mu_k^{(i)}\right), \tag{2.30} \]

\[ f_3^{(i)} = \frac{1}{B_i} \left[ \operatorname{Arcsh}(A_i+B_i) - \operatorname{Arcsh}\left(A_i+B_i\mu_k^{(i)}\right) \right] - \mu_k^{(i)} \cdot f_1^{(i)}, \tag{2.31} \]

\[ f_4^{(i)} = \frac{1}{B_i} \left[ \sqrt{1+(A_i+B_i)^2} - \sqrt{1+\left(A_i+B_i\mu_k^{(i)}\right)^2} \right] - \mu_k^{(i)} \cdot f_2^{(i)} . \tag{2.32} \]

\[ (i=1,2,3,\ldots,n) \]

Let us note, moreover, that \(v_0^{(0)}=0\). The constants \(A_i\) and \(B_i\) entering formulas (2.29)—(2.32), by virtue of the continuity of the pitch-control program, and also by virtue of formulas (2.22)—(2.24), are connected with one another by \(2(n-1)\) relations:

\[ A_i+B_i\mu_k^{(i)}=A_{i+1}+B_{i+1} \tag{2.33} \]

\[ \frac{B_i}{B_{i+1}}=-\frac{v_0^{(i)}}{v_0^{(i+1)}}\cdot\frac{c_i}{c_{i+1}}, \tag{2.34} \]

\[ i=1,2,\ldots,(n-1) \]

Relations (2.34) are obtained in the following way. According to (2.1)

\[ \frac{d}{dt}(\tg\varphi)=-C_2=\mathrm{const}. \]

On the other hand, according to (2.17) and (2.22)

\[ \frac{d}{dt}(\tg\varphi)=\frac{d}{d\mu^{(i)}}(\tg\varphi)\frac{d\mu^{(i)}}{dt} =-B_i\cdot\frac{1-\mu_k^{(i)}}{T_i}=-C_2. \tag{2.35} \]

From (2.35), taking into account (2.21), it is easy to obtain (2.34).

Formulas (2.21)—(2.25) make it possible to calculate the extremal motion completely. In order to obtain a motion satisfying the prescribed isoperimetric conditions, it is necessary to choose appropriately two arbitrary constants of the control program \(C_1\) and \(C_2\), or, what is the same thing, any two constants from the set of \(2n\) quantities \(A_i\) and \(B_i\), for example \(A_1\) and \(B_1\), or \(A_n\) and \(B_n\). Equating \(w\) to zero, and \(u\) to the prescribed value \(u_k\), we obtain, for determining the constants \(A_i\) and \(B_i\), a system of two transcendental equations, whose solution can be found only by trial. It is therefore simpler not to determine the quantities \(A_i\) and \(B_i\), but to prescribe them, and then to calculate all the other quantities from them. Varying \(\mu_k^{(i)}\), \(A_i\), and \(B_i\) within some reasonable range, one can obtain motions with different values of \(v_0^{(i)}\), \(y_k\), and \(u_k\). In this case the computations will be much more economical, since it is not necessary to solve equations, but only to carry out calculations by formulas.

In what follows we shall dwell on the consideration of one very interesting special case of a multistage rocket scheme and shall illustrate on it the method described above.

Let us consider a multistage rocket in which all \(n\) stages are similar to one another in a number of their basic characteristics, namely, a rocket for which

\[ c_i=c;\qquad \mu_k^{(i)}=\mu_k;\qquad v_0^{(i)}=v_0. \tag{2.36} \]

As arbitrary constants of the control program we take \(A_n\) and \(B_n\). All the remaining constants \(A_i\) and \(B_i\) can very simply be expressed in terms of \(A_n\) and \(B_n\). Indeed, according to (2.34)

\[ B_i=B_n, \tag{2.37} \]

and from (2.33), also using (2.37), we easily obtain

\[ A_i=A_n+B_n(1-\mu_k)(n-i)\qquad i=1,2,\ldots,n. \tag{2.38} \]

The calculation scheme will now be as follows. Having prescribed \(\mu_k\), \(A_n\), and \(B_n\) and using formula (2.38), we determine, by formulas (2.29)—(2.32), the quantities \(f_1^{(i)}\), \(f_2^{(i)}\),

\(f_3^{(i)}\) and \(f_4^{(i)}\). From the condition \(u_k=0\) at the end of the motion we have the formula

\[ v_0=\frac{\displaystyle\sum_{i=1}^{n} f_2^{(i)}}{n(1-\mu_k)} . \tag{2,39} \]

On the basis of this formula we determine \(v_0\). After this, by formulas (2,25), (2,27), and (2,28), we determine the quantities \(x_k\), \(y_k\), and \(u_k\).

In carrying out such a calculation, one need not specify the quantities \(A_n\) and \(B_n\) themselves, but rather certain quantities dependent on them.

Fig. 6.

Thus, for example, one may specify the value of the pitch angle at the end of the trajectory and the value of the mass ratio, and compute a series of trajectories, varying the value of the pitch angle at the beginning of the acceleration segment. Denoting the initial pitch angle of the first stage \(\varphi_0^{(1)}\) by \(\varphi_0\), in accordance with (2,23) we shall have

\[ \operatorname{tg}\varphi_0=A_1+B_1 . \tag{2,40} \]

According to (2,24) we obtain

\[ \operatorname{tg}\varphi_k=A_n+B_n\mu_k, \tag{2,41} \]

where \(\varphi_k\) denotes the final pitch angle for the last \(n\)-th stage,

tion \(\varphi_k^{(n)}\). Finally, taking into account that, according to formula (2.38),

\[ A_1=A_n+B_n(1-\mu_k)(n-1) \tag{2.42} \]

and that \(B_1=B_n\), from relations (2.40), (2.41), and (2.42) we obtain formulas relating \(A_n\) and \(B_n\) to \(\varphi_0\) and \(\varphi_k\):

\[ B_n=\frac{\operatorname{tg}\varphi_0-\operatorname{tg}\varphi_k}{n(1-\mu_k)}, \]

\[ A_n=\operatorname{tg}\varphi_k-\mu_k B. \]

The results of the calculation are presented in Figs. 6, 7, and 8. Here the calculations were carried out for two-stage (Fig. 6), three-stage

Fig. 7.

(Fig. 7), and four-stage (Fig. 8) rockets. In the graphs, along the abscissa axis, both in dimensional and dimensionless form (see below), the final velocity of the rocket \(u_k\) is plotted; along the ordinate axis—the final injection altitude \(y_k\). In performing the calculations, the values of the pitch angles at the end of the trajectory were taken equal to

\[ \varphi_k=0^\circ;\ -10^\circ;\ -20^\circ;\ -30^\circ;\ -40^\circ;\ -50^\circ;\ -60^\circ . \]

The values of \(\mu_k\) were taken as

\[ \mu_k = 0.1;\ 0.2;\ 0.3;\ 0.4;\ 0.5. \]

For these parameter values, curves were constructed along which \(\varphi_0\) varied. The extreme right-hand point of each curve corresponds to \(\varphi_0 = 60^\circ\). Exceptions are the curves for which the points corresponding to \(\varphi_0 = 60^\circ\) lie below the abscissa axis; the values of \(\varphi_0\) for the extreme right-hand points of these curves are indicated on the graphs. The remaining points \(\varphi_0\), marked on the graphs by small circles, are plotted at intervals of \(\Delta \varphi_0 = 2^\circ\)

Fig. 8.

in the direction of increasing \(\varphi_0\). It turned out that along such curves the values of \(\nu_0\) change monotonically. For small \(\nu_0\) we obtain points in the lower right-hand part of the curve, corresponding to smaller values of \(\varphi_0\). With increasing \(\nu_0\), we obtain points situated to the left and higher up and corresponding to larger values of the angle \(\varphi_0\). Since, for practically interesting cases, \(\nu_0\) should not be less than unity, the point on the curve

\(\varphi_k = \mathrm{const}\), corresponding to \(\nu_0 = 1\), is the left boundary of the curve under consideration. In Figs. 6–8 the lines connecting the points with \(\nu_0 = 1\), or, equivalently, connecting the left boundaries of the curves \(\varphi_k = \mathrm{const}\), are drawn with a dashed line. It is not difficult to observe that the curves \(\varphi_k = \mathrm{const}\) may lie both to the right and to the left of the curve \(\nu_0 = 1\).

In Figs. 6–8 the curves corresponding to other constant values of \(\nu_0\) are also shown by dashed lines:

\[ \nu_0 = 0.3;\ 0.5;\ 0.7. \]

It is interesting to note that the curves of constant values of \(\nu_0\) may, in a number of cases, intersect. Thus, from Figs. 6–8 it is easy to see that in the cases \(n = 2\), \(\mu_k = 0.3\); \(n = 3\), \(\mu_k = 0.4\); \(n = 4\), \(\mu_k = 0.5\), the curve \(\nu_0 = 1\) intersects the curves \(\nu_0 = 0.7\), \(\nu_0 = 0.5\), and, apparently, the curve \(\nu_0 = 0.3\) (if the curve \(\nu_0 = 1\) is continued in the direction of positive values of \(\varphi_k\)). This means that, for given \(n\) and \(\mu_k\), and for two different values of \(\nu_0\), the optimal programs corresponding to these values of \(\nu_0\) can ensure the attainment of one and the same combination \((y_k u_k)\).

In order to explain the significance of this fact, let us give an example. Suppose that for \(n = 2\), \(\mu_k = 0.3\), \(\nu_0 = 0.5\), and for the altitude \(\tilde y_k = 1.53\) (taken in dimensionless form), it is necessary to find the optimal program ensuring attainment of the maximum velocity \(\tilde u_k\).

Using Fig. 6, by interpolation we find the parameters of the optimal program \(\varphi_0\) and \(\varphi_k\):

\[ \varphi_0 \approx 68^\circ,\quad \varphi_k \approx -43^\circ. \]

The optimal program found ensures attainment of the maximum velocity \(\tilde u_k\), and in this case \(\tilde u_k = 0.46\) (also taken in dimensionless form). However, from the same Fig. 6 it is easy to see that the same combination \((y_k, u_k)\) can be obtained by using another value of \(\nu_0\), namely, \(\nu_0 = 1\). In this case the parameters of the optimal program \(\varphi_0\) and \(\varphi_k\) will already be different,

\[ \varphi_0 = 64^\circ,\quad \varphi_k = 0. \]

It should be borne in mind that, with the initial weight of the rocket kept fixed, larger values of \(\nu_0\), according to formula (2.18), will correspond to smaller values of the rocket-engine thrust \(P\) and, consequently, to a smaller weight of the engine itself. Hence, keeping \(n\) and \(\mu_k\) fixed and ensuring the value \(\nu_0 = 1\) instead of \(\nu_0 = 0.5\), we have the possibility of placing a larger payload in the rocket, without thereby impairing the principal flight characteristics of the rocket—\(y_k\) and \(u_k\).

From the graphs given we also see that, for a given number of rocket stages and for a given \(\mu_k\), there are certain limitations on the possibility of obtaining various combinations of altitude and injection velocity. These limitations have a dual character. To show this, let us divide all the curves \(\varphi_k = \mathrm{const}\) into two classes. To the first class we assign curves having portions situated to the right and above the curve \(\nu_0 = 1\). This class includes, for example, all curves \(\varphi_k = \mathrm{const}\) corresponding to the case \(n = 2\), \(\mu_k = 0.3\), or, for example, the curves \(\varphi_k = -20^\circ\) and \(\varphi_k = -10^\circ\), corresponding to the case \(n = 2\), \(\mu_k = 0.2\) (Fig. 6). To the second class we assign curves located entirely to the left and below the dashed curve \(\nu_0 = 1\). This class includes, for example, all curves \(\varphi_k = \mathrm{const}\) corresponding to the case \(n = 2\), \(\mu_k = 0.1\), or, for example, the curves \(\varphi_k = -60^\circ;\ -50^\circ;\ -40^\circ;\ -30^\circ\), corresponding to the case \(n = 2\), \(\mu_k = 0.2\) (Fig. 6).

From Figs. 6–8 it is easy to see that, for families of curves of the first class, one can easily construct envelopes tangent to curves with identical values of \(\mu_k\) (these envelopes are not shown in the figures). These envelopes bound, in the \(u_k, y_k\) plane, certain regions of attainability; moreover, this limitation is of an energy character and is connected with the possibility, for a given fuel supply, of obtaining only certain limited altitudes and velocities.

Another limitation, pertaining entirely to the curves of the second class, is connected with the practical inexpediency of using rockets for which the thrust of the first stage is less than its initial weight. In this case the dashed curve \(\nu_0=1\) will be the limiting curve in the \((y_k,u_k)\) plane, above which we cannot rise if \(n\) and \(\mu_k\) are given.

Combining the envelope of the family of curves of the first class and the limiting curve of the family of curves of the second class, one can obtain an overall limiting curve indicating the region of attainable combinations of \(y_k\) and \(u_k\) for given \(n\) and \(\mu_k\). The indicated limiting curve gives the greatest values of the velocity attainable at a given altitude for given \(n\) and \(\mu_k\), or the greatest values of the altitude for a given terminal velocity for given \(n\) and \(\mu_k\). In addition, these limiting curves indicate the greatest values of \(\mu_k\) for which the attainment of a given altitude and velocity is energetically possible. Points of the curves \(\mu_k=\mathrm{const}\) lying on the limiting curve will correspond to the energetically most advantageous values of \(\nu_0\).

In Figs. 6–8 the principal scale is the scale for the dimensionless velocity and altitude \(\widetilde{u}_k\) and \(\widetilde{y}_k\), where

\[ \widetilde{u}_k=\frac{u_k}{c},\qquad \widetilde{y}_k=\frac{y_k}{\dfrac{c^2}{g}}. \]

Thanks to this, the graphs presented have a universal form and are suitable for use in the case of various exhaust velocities of the gases from the jet engine nozzle. On the same graphs three dimensional scales are plotted, corresponding to three different exhaust velocities

\[ c=2.5\ \text{km/sec};\quad 3\ \text{km/sec};\quad 3.5\ \text{km/sec}. \]

For the exhaust velocities given, part of the curves falls in the region of such large altitudes \(y_k\) that the direct application of the developed theory of launching in a plane-parallel gravitational field becomes difficult. Nevertheless, even in this case the indicated curves make it possible to draw certain qualitative conclusions of a general character concerning the influence of the choice of rocket parameters and of its pitch-control program on the expansion of the region of attainable velocities and altitudes.

§ 3. THE PROBLEM OF LAUNCHING INTO ORBIT WITH ACCOUNT TAKEN OF THE VARIABILITY OF THE GRAVITATIONAL FIELD AND THE ROTATION OF THE EARTH

In the present section a more complicated formulation of the problem of launching into orbit is considered, taking into account the variability of the gravitational field and the rotation of the Earth. In doing so, the relative motion of the rocket is considered in a coordinate system connected with the Earth.

Assuming that the trajectory of the rocket relative to the Earth is plane, let us consider a Cartesian coordinate system whose origin is placed at the launch point, and whose \(x\)-axis is directed along the tangent to the Earth in the direction

of the rocket’s motion; the \(y\)-axis is directed vertically upward; the \(z\)-axis is perpendicular to the \(x, y\) axes and is directed so as to form with them a right-handed coordinate system (Fig. 9). The equations of motion of the rocket in projections onto the axes of the Cartesian coordinate system are written in the form:

\[ \ddot{x}=-g_x+p_x-a_{\mathrm{per}\,x}-a_{\mathrm{cor}\,x}, \tag{3,1} \]

\[ \ddot{y}=-g_y+p_y-a_{\mathrm{per}\,y}-a_{\mathrm{cor}\,y}, \tag{3,2} \]

\[ \ddot{z}=-g_z+p_z-a_{\mathrm{per}\,z}-a_{\mathrm{cor}\,z}, \tag{3,3} \]

where \(g_x, g_y, g_z\) are the projections of the acceleration of gravity, and \(p_x, p_y, p_z\) are the projections of the acceleration of the reactive force onto the axes \(x, y, z\); \(a_{\mathrm{per}\,x}, a_{\mathrm{per}\,y}, a_{\mathrm{per}\,z}; a_{\mathrm{cor}\,x}, a_{\mathrm{cor}\,y}, a_{\mathrm{cor}\,z}\) are, respectively, the projections of the translational acceleration and the Coriolis acceleration onto the same axes.

Let us analyze the separate terms of equations (3,1)—(3,3). First of all, it is clear that, by virtue of the condition that the rocket move in the plane \(x, y\), one must have

\[ \ddot{z}=\dot{z}=z=0. \tag{3,4} \]

Further, by virtue of the same condition, and also by virtue of the centrality of the gravitational field under consideration, the projection of the acceleration of gravity onto the \(z\)-axis must be equal to zero,

\[ g_z=0. \tag{3,5} \]

Fig. 9.

The projections of the acceleration of gravity onto the axes \(x\) and \(y\) may be represented in the form

\[ g_x=\frac{gR^2}{x^2+(y+R)^2}\cdot \frac{x}{\sqrt{x^2+(y+R)^2}}, \tag{3,6} \]

\[ g_y=\frac{gR^2}{x^2+(y+R)^2}\cdot \frac{y+R}{\sqrt{x^2+(y+R)^2}}. \tag{3,7} \]

The projections of the acceleration of the reactive force are equal to

\[ p_x=p\cos\varphi\cos\beta, \tag{3,8} \]

\[ p_y=p\sin\varphi, \tag{3,9} \]

\[ p_z=p\cos\varphi\sin\beta, \tag{3,10} \]

where \(\varphi\) is the angle between the longitudinal axis of the rocket and the plane \(x, z\) (the pitch angle), \(\beta\) is the angle between the vertical plane passing through the rocket axis and the plane \(x, y\) (the yaw angle), and \(p\), as before, is the magnitude of the acceleration of the reactive force, a prescribed function of time.

The projections of the Coriolis acceleration \(a_{\mathrm{cor}\,x}, a_{\mathrm{cor}\,y}\) and \(a_{\mathrm{cor}\,z}\) are easily expressed in terms of the kinematic parameters of the relative motion if one uses the formula

\[ \mathbf{a}_{\mathrm{cor}}=2[\boldsymbol{\omega}\times \mathbf{v}_{\mathrm{rel}}], \tag{3,11} \]

where \(\boldsymbol{\omega}\) is the vector of the angular velocity of the Earth’s rotation, and \(\mathbf{v}_{\mathrm{rel}}\) is the relative velocity of the rocket’s motion. Note that for the projections of the angular velocity of rotation \(\boldsymbol{\omega}\) onto the axes \(x, y\), and \(z\) we have the formulas

\[ \omega_x=\omega\cos\psi_0\cos\alpha_0, \tag{3,12} \]

\[ \omega_y=\omega\sin\psi_0, \tag{3,13} \]

\[ \omega_z=-\omega\cos\psi_0\sin\alpha_0, \tag{3,14} \]

where \(\omega\) is the magnitude of the angular velocity of the Earth’s rotation, \(\psi_0\) is the latitude of the launch point, and \(\alpha_0\) is the angle between the plane \(x, y\) and the meridional plane.

(azimuth of the trajectory plane). Taking into account (3.11) for the projections of the Coriolis acceleration, we obtain the well-known formulas

\[ a_{\mathrm{cor}\,x}=2[\boldsymbol{\omega}\times \mathbf{v}_{\mathrm{rel}}]_x =\dot{y}\omega_z-z\omega_y, \tag{3.15} \]

\[ a_{\mathrm{cor}\,y}=2[\boldsymbol{\omega}\times \mathbf{v}_{\mathrm{rel}}]_y =z\omega_x-\dot{x}\omega_z, \tag{3.16} \]

\[ a_{\mathrm{cor}\,z}=2[\boldsymbol{\omega}\times \mathbf{v}_{\mathrm{rel}}]_z =\dot{x}\omega_y-\dot{y}\omega_x. \tag{3.17} \]

Finally, taking into account that, by virtue of the conditions of motion, \(z=0\), we shall present the formulas for the Coriolis projections in the following final form:

\[ a_{\mathrm{cor}\,x}=\dot{y}\omega_z, \tag{3.18} \]

\[ a_{\mathrm{cor}\,y}=-\dot{x}\omega_z, \tag{3.19} \]

\[ a_{\mathrm{cor}\,z}=\dot{x}\omega_y-\dot{y}\omega_x. \tag{3.20} \]

Turning now to the projections of the transport acceleration, we note that they are negligibly small in comparison with the other terms of equations (3.1)—(3.3), and in solving the problem under consideration they may be omitted.

Taking into account the remarks made, we write the system of equations (3.1)—(3.3) in the form

\[ \ddot{x}=-g_x+p\cos\varphi\cos\beta-\dot{y}\omega_z, \tag{3.21} \]

\[ \ddot{y}=-g_y+p\sin\varphi+\dot{x}\omega_z, \tag{3.22} \]

\[ 0=p\cos\varphi\sin\beta-(\dot{x}\omega_y-\dot{y}\omega_x). \tag{3.23} \]

The last equation of the system (3.21)—(3.23) serves to determine the angle \(\beta\), necessary for compensating the lateral forces and ensuring the motion of the rocket in the plane \(x,y\). According to (3.23) we have

\[ \sin\beta=\frac{\dot{x}\omega_y-\dot{y}\omega_x}{p\cos\varphi}. \tag{3.24} \]

An estimate of the angle \(\beta\), made with the aid of formula (3.24), shows that in cases of practical interest the angle \(\beta\) will be sufficiently small and, in any case, will never exceed several degrees, so that, in solving the problem under consideration, without great error one may put \(\cos\beta=1\) in equation (3.21).

Taking into account that in the problem under consideration the extent of the ascent trajectory is assumed to be more or less limited, in order to facilitate the solution of the problem we shall simplify the expressions for the components of the acceleration of terrestrial gravity \(g_x\) and \(g_y\), replacing formulas (3.6) and (3.7) by the corresponding approximate formulas. To this end we expand the functions \(g_x(x,y)\) and \(g_y(x,y)\) in the neighborhood of the point \((0,0)\) in a series in powers of \(x\) and \(y\)

\[ g_x=g_x(0,0)+\left(\frac{\partial g_x}{\partial x}\right)_0\cdot x+ \left(\frac{\partial g_y}{\partial y}\right)_0\cdot y+\ldots, \tag{3.25} \]

\[ g_y=g_y(0,0)+\left(\frac{\partial g_y}{\partial x}\right)_0\cdot x+ \left(\frac{\partial g_y}{\partial y}\right)_0\cdot y+\ldots, \tag{3.26} \]

where for the coefficients of the expansion we shall have the following formulas:

\[ g_x(0,0)=0,\qquad \left(\frac{\partial g_x}{\partial x}\right)_0=\frac{g}{R},\qquad \left(\frac{\partial g_x}{\partial y}\right)_0=0, \]

\[ g_y(0,0)=g,\qquad \left(\frac{\partial g_y}{\partial x}\right)_0=0,\qquad \left(\frac{\partial g_y}{\partial y}\right)_0=-\frac{2g}{R}. \]

In what follows, in the expansions (3.25), (3.26) we shall confine ourselves to terms containing \(x\) and \(y\) to the first degree. We obtain the following approximate

formulas:

\[ g_x = g \cdot \frac{x}{R}, \tag{3.27} \]

\[ g_y = g - 2g \frac{y}{2}. \tag{3.28} \]

Taking into account all the remarks and simplifications made, we shall present the system of equations of motion of the rocket in the following final form:

\[ \psi_1 = \dot u - 2\omega_z w + \nu^2 x - p \cos \varphi = 0, \tag{3.29} \]

\[ \psi_2 = \dot w + 2\omega_z u - 2\nu^2 y - p \cos \varphi + g = 0, \tag{3.30} \]

\[ \psi_3 = \dot x - u = 0, \tag{3.31} \]

\[ \psi_4 = \dot y - w = 0, \tag{3.32} \]

where

\[ \nu = \sqrt{\frac{g}{R}}, \tag{3.33} \]

\[ \omega_z = -\omega \cos \psi_0 \sin \alpha_0. \tag{3.34} \]

As in § 1, let us pose the variational problem of attaining the maximum velocity at a prescribed altitude under the condition that the angle between the vector of the final velocity and the local horizon is zero. As before, by virtue of reciprocity the solution obtained will ensure the attainment of the maximum altitude at a prescribed velocity, as well as the attainment of prescribed values of altitude and velocity with minimum fuel consumption.

Let us formulate the boundary conditions. Suppose that at the beginning of the motion (at \(t=t_0\)) we have the coordinates \(x_0\) and \(y_0\) and certain values of the horizontal and vertical projections of the velocity, \(u_0\) and \(w_0\). Let us also assume that the kinematic parameters \(x_0, y_0, u_0\), and \(w_0\) depend on a certain parameter, which may be, in particular, the angle \(\theta_0\) between the vector of the initial velocity \(v_0\) and the horizon at the instant \(t=t_0\). Introducing in the subsequent formulas the angle \(\theta_0\) as such a parameter, we shall keep in mind the possibility of replacing \(\theta_0\) by some other parameter. We assume that at \(t=t_0\)

\[ x_0 = x_0(\theta_0), \quad y_0 = y_0(\theta_0), \quad u_0 = u_0(\theta_0), \quad w_0 = w_0(\theta_0). \tag{3.35} \]

At the end of the motion, at the instant \(t=T\), the total flight altitude \(h\) is fixed, and the angle between the velocity vector and the local horizon must be zero. Both of these conditions are written in the following form: for \(t=T\)

\[ \sqrt{(R+y_k)^2 + x_k^2} - R = h, \tag{3.36} \]

\[ \frac{x_k}{R+y_k} = -\frac{w_k}{u_k}. \tag{3.37} \]

In the limiting case, if \(x_k\) is small, condition (3.37) takes the form

\[ w_k = 0, \]

i.e., the same as in the problem considered in §§ 1 and 2.

Thus, let us pose the variational problem. For a prescribed function \(p(t)\), find the function \(\varphi(t)\) that provides the maximum horizontal velocity \(v_k\) at the prescribed altitude \(h\), as well as the most advantageous angle \(\theta_0\) between the vector of the initial velocity and the horizon, corresponding to the point of launch. We note that the dependence of the functional of the problem \(v_k\) on the function \(\varphi(t)\) cannot, as it was earlier, be represented in the simplest form.

D. E. OKHOTSIMSKII, T. M. ENEEV

To solve the problem posed, it is necessary, first of all, to find the first variation of the functional \(v_k\). To find the first variation we shall use the method of Lagrange multipliers. For this purpose, along with the original functional

\[ v_k=\sqrt{u_k^2+w_k^2} \tag{3.38} \]

we consider the new functional

\[ J=\sqrt{u_k^2+w_k^2}+\int_0^T H\,dt. \tag{3.39} \]

For the integrand \(H\) we have the expression

\[ H=\lambda_1\psi_1+\lambda_2\psi_2+\lambda_3\psi_3+\lambda_4\psi_4, \tag{3.40} \]

where \(\lambda_1,\lambda_2,\lambda_3\) and \(\lambda_4\) are certain, as yet undetermined, functions of time, while \(\psi_1,\psi_2,\psi_3\) and \(\psi_4\) are determined according to (3.29)—(3.32). It is clear that the functionals (3.39) and (3.38) coincide by virtue of the fulfillment of conditions (3.29)—(3.32). The variations of the functionals (3.38) and (3.39) will also coincide if, when the function \(\varphi\) is varied, the constraints (3.29)—(3.32) are satisfied, i.e., the conditions

\[ \delta\psi_1=\delta\psi_2=\delta\psi_3=\delta\psi_4=0. \]

are fulfilled. The first variation of the functional \(J\) will have the form

\[ \delta J= \frac{u_k}{\sqrt{u_k^2+w_k^2}}\,\delta u_k+ \frac{w_k}{\sqrt{u_k^2+w_k^2}}\,\delta w_k+ \int_0^T \delta H\,dt, \tag{3.41} \]

where

\[ \delta H=H'_u\delta u+H'_{\dot u}\delta\dot u+ H'_w\delta w+H'_{\dot w}\delta\dot w+ H'_x\delta x+H'_{\dot x}\delta\dot x+ \]

\[ +H'_y\delta y+H'_{\dot y}\delta\dot y+H'_\varphi\delta\varphi. \tag{3.42} \]

For the derivatives \(H'_u\), \(H'_{\dot u}\), etc., according to (3.29)—(3.32) and (3.40), we have the formulas:

\[ \left. \begin{aligned} H'_u&=2\omega_z\lambda_2-\lambda_3; & H'_w&=-2\omega_z\lambda_1-\lambda_4; & H'_x&=v^2\lambda_1; & H'_y&=-2v^2\lambda_2;\\ H'_{\dot u}&=\lambda_1; & H'_{\dot w}&=\lambda_2; & H'_{\dot x}&=\lambda_3; & H'_{\dot y}&=\lambda_4;\\ H'_\varphi&=p(\lambda_1\sin\varphi-\lambda_2\cos\varphi). \end{aligned} \right\} \tag{3.43} \]

Assuming the functions \(\lambda_1,\lambda_2,\lambda_3\) and \(\lambda_4\) to be continuous and differentiable at every point, we integrate by parts in the integrand expression on the right-hand side of (3.41) the terms containing derivatives of the variations. As a result we obtain

\[ \delta J= \frac{u_k}{v_k}\delta u_k+ \frac{w_k}{v_k}\delta w_k+ \lambda_{1k}\delta u_k+ \lambda_{2k}\delta w_k+ \lambda_{3k}\delta x_k+ \lambda_{4k}\delta y_k- \]

\[ -\lambda_{10}\delta u_0-\lambda_{20}\delta w_0 -\lambda_{30}\delta x_0-\lambda_{40}\delta y_0+ \int_0^T \widetilde{\delta H}\,dt, \tag{3.44} \]

where \(\lambda_{10},\lambda_{20},\lambda_{30},\lambda_{40}\) and \(\lambda_{1k},\lambda_{2k},\lambda_{3k},\lambda_{4k}\) are the values of the functions \(\lambda_1,\lambda_2,\lambda_3\) and \(\lambda_4\), respectively, at the initial and final instants of time. We note that the variations \(\delta u_k,\delta w_k,\delta x_k\) and \(\delta y_k\), by virtue of conditions (3.36) and (3.37), are not independent, but are connected with one another by the relations

\[ \frac{1}{R+y_k}\delta x_k- \frac{x_k}{(R+y_k)^2}\delta y_k = \frac{w_k}{u_k^2}\delta u_k-\frac{1}{u_k}\delta w_k, \tag{3.45} \]

\[ x_k\delta x_k+(R+y_k)\delta y_k=0. \tag{3.46} \]

Thus, of the four variations \(\delta u_k,\delta w_k,\delta x_k\) and \(\delta y_k\), only two are independent. In what follows, as independent variations

we take the variations \(\delta u_k\) and \(\delta w_k\). Also, according to (3.35), the variations \(\delta u_0,\ \delta w_0,\ \delta x_0\) and \(\delta y_0\) are not independent, but are related to one another through the variation \(\delta\theta_0\).

Obviously, we have

\[ \delta u_0=\frac{\partial u_0}{\partial \theta_0}\delta\theta_0,\quad \delta w_0=\frac{\partial w_0}{\partial \theta_0}\delta\theta_0,\quad \delta x_0=\frac{\partial x_0}{\partial \theta_0}\delta\theta_0,\quad \delta y_0=\frac{\partial y_0}{\partial \theta_0}\delta\theta_0. \tag{3.47} \]

Using (3.45), (3.46), and (3.47), let us eliminate from the expression preceding the integral in the right-hand side of (3.44) the superfluous variations. After simple transformations we obtain:

\[ \begin{aligned} \delta J={}& \left\{ \frac{u_k}{v_k}+\lambda_{1k} +\frac{w_k}{u_k^2}\left(\frac{R+y_k}{R+h}\right)^2 \left[(R+y_k)\lambda_{3k}+x_k\lambda_{4k}\right] \right\}\delta u_k+ \\ &+ \left\{ \frac{w_k}{v_k}+\lambda_{2k} -\frac{1}{u_k}\left(\frac{R+y_k}{R+h}\right)^2 \left[(R+y_k)\lambda_{3k}+x_k\lambda_{4k}\right] \right\}\delta w_k+ \\ &+ \left\{ \lambda_{10}\frac{\partial u_0}{\partial\theta_0} +\lambda_{20}\frac{\partial w_0}{\partial\theta_0} +\lambda_{30}\frac{\partial x_0}{\partial\theta_0} +\lambda_{40}\frac{\partial y_0}{\partial\theta_0} \right\}\delta\theta_0 +\int_0^T \widetilde H\,dt. \end{aligned} \tag{3.48} \]

Further, for the expression \(\delta\widetilde H\) standing under the integral in (3.48), we have the formula

\[ \begin{aligned} \delta\widetilde H={}& (2\omega_z\lambda_2-\lambda_3-\dot\lambda_1)\delta u +(-2\omega_z\lambda_1-\lambda_4-\dot\lambda_2)\delta w \\ &+(\nu^2\lambda_1-\dot\lambda_3)\delta x +(-2\nu^2\lambda_2-\dot\lambda_4)\delta y +p(\lambda_1\sin\varphi-\lambda_2\cos\varphi)\delta\varphi . \end{aligned} \tag{3.49} \]

Let us now choose the multipliers \(\lambda_1,\lambda_2,\lambda_3\), and \(\lambda_4\) in such a way that, in the right-hand side of the expression for the variation of the functional (3.48), the terms under the integral containing the variations \(\delta u,\delta w,\delta x,\delta y\), as well as the terms standing before the integral and containing the variations \(\delta u_k\) and \(\delta w_k\), vanish. As a result we obtain the system of differential equations for \(\lambda_1,\lambda_2,\lambda_3\), and \(\lambda_4\):

\[ \dot\lambda_1=2\omega_z\lambda_2-\lambda_3, \tag{3.50} \]

\[ \dot\lambda_2=-2\omega_z\lambda_1-\lambda_4, \tag{3.51} \]

\[ \dot\lambda_3=\nu^2\lambda_1, \tag{3.52} \]

\[ \dot\lambda_4=-2\nu^2\lambda_2, \tag{3.53} \]

and also two conditions at the right boundary, which we shall present in the form

\[ \lambda_{1k}=-\frac{u_k}{v_k} -\frac{w_k}{u_k^2}\left(\frac{R+y_k}{R+h}\right)^2 \left[(R+y_k)\lambda_{3k}+x_k\lambda_{4k}\right], \tag{3.54} \]

\[ \lambda_{2k}=-\frac{w_k}{v_k} +\frac{1}{u_k}\left(\frac{R+y_k}{R+h}\right)^2 \left[(R+y_k)\lambda_{3k}+x_k\lambda_{4k}\right]. \tag{3.55} \]

The first variation of the functional (3.39) can now be represented in a form that directly contains the connection of the variation of the functional with the variations of the angle \(\theta_0\) and of the function \(\varphi(t)\):

\[ \delta J= \left( \lambda_{10}\frac{\partial u_0}{\partial\theta_0} +\lambda_{20}\frac{\partial w_0}{\partial\theta_0} +\lambda_{30}\frac{\partial x_0}{\partial\theta_0} +\lambda_{40}\frac{\partial y_0}{\partial\theta_0} \right)\delta\theta_0 +\int_0^T p(\lambda_1\sin\varphi-\lambda_2\cos\varphi)\delta\varphi\,dt. \tag{3.56} \]

Assuming that there exists an interior extremum of the functional (3.39) both with respect to \(\theta_0\) and with respect to \(\varphi(t)\), we find it from the condition

\[ \delta J=0. \tag{3.57} \]

Since the variations \(\delta\theta_0\) and \(\delta\varphi\) are independent, condition (3.57) will obviously be satisfied if

\[ \lambda_{10}\frac{\partial u_0}{\partial\theta_0} +\lambda_{20}\frac{\partial w_0}{\partial\theta_0} +\lambda_{30}\frac{\partial x_0}{\partial\theta_0} +\lambda_{40}\frac{\partial y_0}{\partial\theta_0}=0, \tag{3.58} \]

\[ \tg\varphi=\frac{\lambda_2}{\lambda_1}. \tag{3.59} \]

D. E. OKHOTSIMSKY, T. M. ENEEV

We have obtained two relations, of which the first serves to determine the optimal angle \(\theta_0\), and the second to determine the optimal program \(\varphi(t)\).

The system of equations (3.50)—(3.53), (3.59), together with the boundary conditions (3.54), (3.55), (3.58), and also with the equations of motion (3.29)—(3.32) and the boundary conditions (3.35)—(3.37), gives the complete solution of the problem posed. In this case the dependence of the optimal control program on time can be obtained in explicit form.

Indeed, the system (3.50)—(3.53) is a system of linear differential equations with constant coefficients and therefore can be easily integrated. To find the general solution of the system (3.50)—(3.53) it is necessary to find the roots of the characteristic equation of this system:

\[ k^4-\left(\nu^2-4\omega_z^2\right)k^2-2\nu^4=0. \tag{3.60} \]

Equation (3.60) will have two real and two imaginary roots

\[ k_{1,2}=\pm \nu_1,\qquad k_{3,4}=\pm i\nu_2, \tag{3.61} \]

where

\[ \nu_1=\frac{\sqrt{2}}{2}\sqrt{\sqrt{\left(\nu^2-4\omega_z^2\right)^2+8\nu^4}+\left(\nu^2-4\omega_z^2\right)}, \tag{3.62} \]

\[ \nu_2=\frac{\sqrt{2}}{2}\sqrt{\sqrt{\left(\nu^2-4\omega_z^2\right)^2+8\nu^4}-\left(\nu^2-4\omega_z^2\right)}. \tag{3.63} \]

Finding the general solution of the system (3.50)—(3.53) and substituting it into (3.59), we obtain the general formula for the optimal pitch-control program, containing the explicit dependence of \(\varphi\) on time \(t\):

\[ \tan\varphi= \frac{ \sigma_3\cosh \nu_1(T-t)-\sigma_4\sinh k_1(T-t)+s_3\cos k_2(T-t)-s_4\sin k_2(T-t) }{ s_1\cosh \nu_1(T-t)-s_3\sinh k_1(T-t)+\sigma_1\cos k_2(T-t)-\sigma_3\sin k_2(T-t) }, \tag{3.64} \]

where

\[ \left. \begin{aligned} \sigma_1&=(\nu_1^2+\nu^2+4\omega_z^2)\lambda_{1k}+2\omega_z\lambda_{4k},\\ \sigma_2&=(\nu_2^2+2\nu^2-4\omega_z^2)\lambda_{2k}+2\omega_z\lambda_{3k},\\ \sigma_3&=\frac{1}{\nu_2}\left[2\omega_z(\nu_1^2-\nu^2+4\omega_z^2)\lambda_{2k} -(\nu_1^2+\nu^2+4\omega_z^2)\lambda_{3k}\right],\\ \sigma_4&=-\frac{1}{\nu_1}\left[2\omega_z(\nu_1^2+\nu^2-4\omega_z^2)\lambda_{1k} +(\nu_1^2+2\nu^2-4\omega_z^2)\lambda_{4k}\right], \end{aligned} \right\} \tag{3.65} \]

\[ \left. \begin{aligned} s_1&=(\nu_2^2-\nu^2-4\omega_z^2)\lambda_{1k}-2\omega_z\lambda_{4k},\\ s_2&=(\nu_1^2-2\nu^2+4\omega_z^2)\lambda_{2k}-2\omega_z\lambda_{3k},\\ s_3&=\frac{1}{\nu_1}\left[2\omega_z(\nu_2^2+\nu^2-4\omega_z^2)\lambda_{2k} -(\nu_2^2-\nu^2-4\omega_z^2)\lambda_{3k}\right],\\ s_4&=-\frac{1}{\nu_2}\left[2\omega_z(\nu_1^2-\nu^2+4\omega_z^2)\lambda_{1k} +(\nu_1^2-2\nu^2+4\omega_z^2)\lambda_{4k}\right]. \end{aligned} \right\} \tag{3.66} \]

In order to carry the problem through to completion, it is necessary in formulas (3.65) and (3.66) to determine the constants \(\lambda_{1k}, \lambda_{2k}, \lambda_{3k}\), and \(\lambda_{4k}\). Let us recall here that \(\lambda_{1k}\) and \(\lambda_{2k}\) are determined through \(\lambda_{3k}, \lambda_{4k}\) and through the kinematic parameters \(u_k, w_k, x_k, y_k\) by formulas (3.54), (3.55).

Thus formula (3.64) determines the control program \(\varphi(t)\) as a function of time, and also of six parameters—\(\lambda_{3k}, \lambda_{4k}, u_k, w_k, x_k, y_k\), i.e.

\[ \varphi=\varphi(t;\lambda_{3k},\lambda_{4k},u_k,w_k,x_k,y_k). \tag{3.67} \]

Substitute \(\varphi(t)\), defined by formula (3.64), into the equations of motion (3.29)—(3.32). The system (3.29)—(3.32) will then be transformed into a nonhomogeneous system of linear differential equations with constant coefficients and with a variable right-hand side, explicitly dependent on time. The solution of such a system can be represented in quadratures. As a result of integrating (3.29)—(3.32) we obtain four relations connecting the seven parameters \(u_k, w_k, x_k, y_k, \lambda_{3k}, \lambda_{4k}\), and \(\theta_0\). We write these relations in the form

\[ u_k-u_0(\theta_0)+\Phi_1(\lambda_{3k},\lambda_{4k},u_k,w_k,x_k,y_k)=0, \tag{3.68} \]

\[ w_k-w_0(\theta_0)+\Phi_2(\lambda_{3k},\lambda_{4k},u_k,w_k,x_k,y_k)=0, \tag{3.69} \]

\[ x_k-x_0(\theta_0)+\Phi_3(\lambda_{3k},\lambda_{4k},u_k,w_k,x_k,y_k)=0, \tag{3.70} \]

\[ y_k-y_0(\theta_0)+\Phi_4(\lambda_{3k},\lambda_{4k},u_k,w_k,x_k,y_k)=0, \tag{3.71} \]

where \(\Phi_1,\Phi_2,\Phi_3\), and \(\Phi_4\) contain integrals which, in the general case, are not expressible in elementary functions. Equations (3.68)—(3.71), together with the terminal relations (3.36)—(3.37) and relation (3.58), constitute a closed system of seven transcendental equations with seven unknowns—\(u_k,w_k,x_k,y_k,\lambda_{3k},\lambda_{4k},\theta_0\)*). Solving this system, we find the values of the kinematic parameters at the end of the injection trajectory \(u_k,w_k,x_k,y_k\), the optimal value of the initial angle \(\theta_0\), and the values of the parameters of the optimal control program \(\lambda_{1k}, \lambda_{2k}, \lambda_{3k}, \lambda_{4k}\); moreover, \(\lambda_{1k}\) and \(\lambda_{2k}\) are computed additionally with the aid of formulas (3.54) and (3.55). The system of equations (3.68)—(3.71), (3.36), (3.37), and (3.58) is very complicated and cumbersome, and in the general case it is impossible to solve it explicitly. Therefore, in practical use of the proposed method for finding the optimal program, it is expedient to employ some numerical iterative method to solve the indicated system. Without dwelling here on the question of practical use of the optimal program found, let us analyze formula (3.64) and simplify it.

First of all, let us clarify the influence of the Earth’s rotation on the choice of the optimal pitch program and show that this influence is insignificant. Consider the case when the conditions of the rocket’s motion are such that

\[ \omega_y=0. \tag{3.72} \]

Condition (3.72) is satisfied, for example, if \(\alpha_0=0\) or \(\pi\), and also if

\[ \psi_0=\pm\frac{\pi}{2}. \]

In this case formulas (3.62) and (3.63) for \(\nu_1\) and \(\nu_2\) are greatly simplified and take the form

\[ \nu_1=\sqrt{2}\,\nu,\qquad \nu_2=\nu. \]

The formulas for \(\sigma_1,\sigma_2,\sigma_3\), and \(\sigma_4\) are also simplified and take the form

\[ \left. \begin{aligned} \sigma_1&=3\nu^2\lambda_{1k}, & \sigma_3&=-3\lambda_{3k},\\ \sigma_2&=3\nu^2\lambda_{2k}, & \sigma_4&=-3\lambda_{4k}\frac{1}{\sqrt{2}}. \end{aligned} \right\} \tag{3.73} \]

As for \(S_1,S_2,S_3\), and \(S_4\), in this case it is easy to obtain from (3.66) that

\[ S_1=S_2=S_3=S_4=0. \tag{3.74} \]

Taking (3.73) and (3.74) into account, we obtain the considerably simplified formula for the program \(\varphi(t)\)

\[ \operatorname{tg}\varphi = \frac{ \lambda_{3k}\operatorname{ch}\sqrt{2}\nu( T-t) + \dfrac{\lambda_{4k}}{\nu\sqrt{2}}\operatorname{sh}\sqrt{2}\nu(T-t) }{ \lambda_{1k}\cos\nu(T-t) + \dfrac{\lambda_{3k}}{\nu}\sin\nu(T-t) }. \tag{3.75} \]

\[ \text{*) Note that } \lambda_{10},\lambda_{20},\lambda_{30},\text{ and }\lambda_{40},\text{ which enter (3.58), are expressed by means of terminal formulas through } \lambda_{3k},\lambda_{4k},u_k,w_k,x_k,\text{ and }y_k. \]

Formula (3.75) can be used, with a high degree of accuracy, as the formula for the optimal program also in the case when \(\omega_z \ne 0\). This is easy to prove if one takes into account that the rotation of the Earth affects the form of the optimal program through the quantities \(\nu_1\) and \(\nu_2\).

From (3.62) and (3.63) we see that \(\omega_z\) enters the formulas for \(\nu_1\) and \(\nu_2\) only in the combination \((\nu^2 - 4\omega_z^2)\). It can be shown that \(4\omega_z^2\) is small in comparison with \(\nu^2\). Using this, with the aid of (3.62) and (3.63), we derive approximate formulas for \(\nu_1\) and \(\nu_2\):

\[ \nu_1 \simeq \sqrt{2}\,\nu \left[1 - \frac{2}{3}\left(\frac{\omega_z}{\nu}\right)^2\right], \qquad \nu_2 \simeq \nu \left[1 + \frac{2}{3}\left(\frac{\omega_z}{\nu}\right)^2\right]. \tag{3.76} \]

Let us take for \(\omega_z\) the maximum possible value, equal to \(\omega\). It is known that

\[ \omega = 7.292 \cdot 10^{-5}. \tag{3.77} \]

On the other hand, by formula (3.33) we compute \(\nu\)

\[ \nu = 1.241 \cdot 10^{-3}. \tag{3.78} \]

Using (3.77), (3.78), we compute the maximum value of the correction term in formulas (3.76):

\[ \frac{2}{3}\left(\frac{\omega}{\nu}\right)^2 = 0.230 \cdot 10^{-2}. \tag{3.79} \]

Thus, the correction to \(\nu_1\) and \(\nu_2\) due to the rotation of the Earth is indeed very small and, taking into account that the variational problem under consideration has been solved with a number of simplifications in the formulas for the external acting forces, this correction may certainly be neglected.

Thus, formula (3.75) may be regarded as a general formula for the optimal pitch-control program, suitable for use under various geographic conditions of injection into orbit.

Let us consider the limiting form of formula (3.75), corresponding to a small value of the parameter \(\nu T\). In this case, putting \(\operatorname{ch}\sqrt{2}\nu T\) and \(\cos \nu T\) equal to unity, and \(\operatorname{sh}\sqrt{2}\nu T\) and \(\sin \nu T\) equal to \(\nu T\), according to (3.75) we shall obviously have

\[ \operatorname{tg}\varphi = \frac{\lambda_{2k} + \lambda_{4k}(T - t)} {\lambda_{1k} + \lambda_{3k}(T - t)}. \tag{3.80} \]

Formula (3.80) can also be obtained if the problem under consideration is solved without taking into account the variability of the gravitational-force field. If, in doing so, the quantities

\[ \frac{x_k}{R}, \qquad \frac{y_k}{R}, \qquad \frac{T}{\sqrt{\dfrac{R}{g}}} \quad \text{and} \quad \frac{w_k}{u_k}, \]

are regarded as sufficiently small, then as a result we obtain formula (2.1) of the preceding section, and together with it the entire theory of optimal injection into orbit developed for a plane-parallel gravitational field.