用于最优控制的简单软件

在本文中，我想介绍一种紧凑且易于使用的软件工具，以实现对动态系统的最佳控制。最佳控制是一个非常广泛的主题。其基础是由Richard Bellman，Lev Pontryagin，Rudolf Kalman等人在20世纪50年代末至60年代奠定的。最佳控制具有许多实际应用。在这里，我们将重点关注以下经典的最优控制问题。我们的目标是通过外部控制输入将给定的动态系统从初始状态转移到某个最终状态。我们希望在此过渡期间系统参数尽可能接近其规定的轨迹。这是具有许多实际应用的最佳控制的非常基本的任务。例如，启动电动机将电流的初始跳跃限制在可接受的值。

处理这个问题甚至它的陈述需要大量的数学。我将仅限于描述问题所需的最低要求，并了解如何使用软件。附录中给出了一些额外的递归公式，没有证明。

问题陈述

让我们正式化我们的问题。动态系统用以下微分方程组描述：

用于最优控制的简单软件

（1）

其中t是时间，x是状态向量，由表征系统行为的参数组成，ẋ是x的一阶导数的向量，m是控制向量。通常，控制矢量的维度不超过状态矢量的维度。在该语句中，向量m独立于向量x，并且可以被认为是开环系统的控制。具有输入控制向量m和状态向量x的开环控制系统如下所示：

用于最优控制的简单软件

开环控制系统

我们希望在0到t f的最终时间之间的时间间隔内获得控制策略m（t），这最小化了成本函数。

用于最优控制的简单软件

（2）

其中上标T表示矢量或矩阵的转置，r是期望的状态向量，u是期望的控制向量，Q和Z分别是用于状态和控制向量的期望参数和实际参数之间的偏差的权重矩阵。在大多数情况下Q和Z.是对角矩阵（但不一定！）。它们的每个对角线元素定义在给定时刻参数的期望值和实际值之间的差的平方的相对权重。矩阵可以是时间依赖的，例如，通常的做法是在最后时间增加差异的相对权重以确保期望的最终状态。这里提出的问题通常被称为“跟踪问题”。

线性动态系统

在我们找到解决方案之前，让我们讨论动态系统的一个特例。研究最充分的系统类是线性动态系统，即用一组线性微分方程描述动力学的系统。对于这样的系统，等式（1）可以重写如下：

用于最优控制的简单软件

（3）

其中A(t) 和B(t) 是矩阵。首先，我们将为这类系统提供解决方案，然后将其扩展到非线性动态系统。

离散时间

由于我们将提供数字解决方案，我们将在离散的时刻处理系统。因此，我们将重新设计离散时间t = k·Δt的问题，其中Δt构成一个小的时间间隔（采样间隔），k是给定的时间步长，0 <= k <N，t f =（N-1）⋅ Δt。现在可以在离散时间内重写等式（1）：

用于最优控制的简单软件

（4）

和线性系统定义为

用于最优控制的简单软件（5）

其中˚F ķ和ħ ķ是从 A(t)和B(t) 中获得的矩阵，如下：
Fk = exp(A(k⋅Δt)⋅Δt) ≅ I+A(k⋅Δt)⋅Δt, Hk = A-1[exp(A(k⋅ Δt)⋅Δt) - I]B ≅ B(k⋅Δt)⋅Δt，其中I是单位矩阵。

离散时间内的成本函数（2）将被重新整形为

用于最优控制的简单软件

（6）

线性案例的解决方案

可以基于Bellman的最优性原理获得问题的数值解。结果非常笨重，并在附录中给出。为了一般理解，必须知道m k控制策略是在两个主要运行中计算的。首先，从最终时间步骤N-1开始并且下降到初始时间0来执行逆计算运行。使用系统矩阵F和H运行逆计算，期望状态r和控制u向量和权重矩阵Q，Z产生一些向量c k和矩阵L k。然后从时间0到N-1执行直接运行，在每个k时间步长获得矢量m k

用于最优控制的简单软件

（7）

反向和直接计算运行如下图所示：

用于最优控制的简单软件

计算方案

向量c N-k和矩阵L N-k从逆运行中已知并且取决于系统矩阵F和H，用户定义的期望向量和权重。它们还依赖于参与反向运行并从零初始条件开始的辅助矩阵P N-k和向量v N-k。在直接计算运行的每个步骤中，我们基于此时的状态向量x k获得控制向量m k。然后使用等式（5）计算下一状态x k + 1。重复该过程直到最后一刻N-1给我们提供了开环控制策略m k和系统x k与时间的状态。声明（7）为我们提供了闭环控制，如下所示

用于最优控制的简单软件

闭环控制系统

其中向量c构成闭环控制，矩阵L提供来自状态向量x的反馈。这些参数允许设计者根据系统的实际状态生成矢量m，从而使系统更加稳定和健壮。

线性案例解决方案的探讨

对于线性系统，控制策略m k最小成本函数（6）仅通过一次迭代包括一次逆运行和一次直接运行来实现。如上所述，该策略取决于系统矩阵和用户选择的参数，即期望的向量和权重矩阵。这些参数足够直观，可以将控制策略设计减少到与它们一起玩。用户选择一组参数并计算控制策略和适当的状态轨迹。如果她/他对结果不满意，则对另一组期望的矢量和权重矩阵重复计算。观察到系统矩阵F，H和权重矩阵Q也是有趣的时间不变，然后反馈矩阵L也将是时不变的，这大大简化了系统控制器的设计。对于线性情况系统矩阵F，H通常是时不变的。因此，为了获得时间不变反馈，用户应选择恒定权重矩阵Q，即在最终时刻没有正式严格实现期望的最终状态的控制。在现实生活中适当选择Q.即使在最后一次之前，我们也能在可接受的时间内达到理想的最终状态。尽管使用状态向量坐标实现控制策略很方便，但在现实生活中通常并非所有这些坐标都可用于测量。已经开发了各种技术来克服这个困难。可以使用不同种类的状态坐标观察器来使用可用于测量的坐标的信息来估计未测量的坐标。

非线性系统

非线性系统案例要复杂得多。与线性系统不同，它不能对各种系统做出总体决策。在这里，我们可以尝试使用拟线性化方法将非线性优化问题简化为线性优化问题。为了使用相同的数学解，我们必须以类似于等式（5）的线性化形式呈现由等式（1）给出的一般动态系统描述。为此，我们应该考虑迭代计算过程。我们可线性（1）至（5）在上的特定点X ķ和Mķ的我个轨迹：

用于最优控制的简单软件

（8）

其中

用于最优控制的简单软件

（9）

线性化系统矩阵F和H是在第i个轨迹上的点k处的函数f的适当偏导数的矩阵。第（i + 1）条轨迹上的状态和控制矢量可以表示为

用于最优控制的简单软件

（10）

可以针对增量Δx和Δm公式化等式（8）。在这种情况下，在成本函数（6）中，新向量r和u分别对应于

用于最优控制的简单软件

（11）

这允许我们将非线性问题减少到一系列拟线性迭代，其中每个都在上面的计算方案图中描述。该算法包括以下步骤：

假设迭代i = 0向量x k和m k（通常所有零接受x k = 0初始条件，但选择取决于用户）根据（8 ）计算线性化的所有F k i和H k i为第i个轨迹）。
根据（11）次迭代为i-th计算新ř ķ i和ü ķ i。
使用反向和正向运行作为线性的情况计算Δm ķ i和ΔX ķ i。
对于第（i + 1）次迭代计算x k i + 1和m k i + 1。
检查条件以停止迭代。该条件可以例如基于与先前迭代相比的成本函数的值之间的差异，或者仅包含固定数量的迭代。无论如何，如果不满足条件，则应重复步骤1-5，直到满足停止条件。在满足条件之后，可以将策略m k和状态x k视为最终解。

获得函数f的偏导数的精确公式可能是繁琐的工作，并且当函数f在数值上定义时变得不可能，例如作为表格。在这些情况下，可以通过给出偏导数的状态向量坐标的小增量来使用简单的数值近似，并计算f函数的结果值。代码示例提供了适当的示例。

方法限制

像任何数字方法一样，这个有其局限性。一个限制是相当明显的，应考虑线性和非线性情况：时间采样间隔应足够小，以保证计算的收敛。r，u和权重因子Q和Z的理想值应该是合理的，以实现直观的最佳系统行为。在非线性系统的情况下，我们应该清楚地理解这是“没有严格规则的游戏”，因此不能保证该技术适用于特定情况。该方法对线性化系统矩阵精度的计算很敏感。可以使用更复杂的方法来估计该方法对特定非线性系统的适用性，但是对于工程师来说，使用尝试失败方法似乎更实用。同样，在某种程度上处理非人体系仍然是“艺术而不是科学”。

不断反馈

与线性情况一样，对于非线性系统，可以考虑具有恒定反馈的控制。为了在非线性情况权重矩阵Q，Z中获得恒定反馈，状态x和控制m的初始轨迹应该是时不变的（初始状态x 0除外），并且应该在一次迭代中实现可接受的结果（因为系统矩阵F和H）根据（9）计算，取决于先前的迭代轨迹）。通常工程师采用相当非正式的方法进行系统优化，当向成本函数提供最小值时，不是最终目标，而是达到可接受系统行为的手段。因此，为了简化反馈，在某些情况下牺牲进一步的成本函数最小化可能是明智的。当然，像非线性系统中的所有内容一样，并非在所有情况下都是如此。下面将讨论适当的数字示例。

软件的使用

该软件使用普通的C＃编写，无需任何其他库。项目MatrixCalc提供矢量和矩阵的算术运算，包括换位和反转。项目SequentialQuadratic支持系统优化和模拟的计算。项目Test使用SequentialQuadratic的静态方法和接口来解决特定问题。DLL SequentialQuadratic有几种内部类型负责不同的优化问题（线性与恒定权重和期望值，线性与时间权重和期望值以及适当的非线性变化）。这些类型的实例是用重载的静态方法SQ.Create()创建的，并通过相应的接口集IBase，IBaseLQExt，ITrackingLQ和INonLinearQ向调用方公开。为了执行优化程序，用户应调用适当的方法提供输入并获得所需的接口。然后应该调用RunOptimization()方法来实际执行优化计算，然后是以逗号分隔格式将结果写入与Excel兼容的* .csv文件的Output2Csv()方法。用给定的控制策略方法Simulate()来模拟系统被调用。用于优化非线性系统（Van der Pol振荡器）的代码示例如下所示：

public static void VanDerPol()
{
	const string TITLE = "Van der Pol";
	Console.WriteLine($"Begin \"{TITLE}\"");

	var lstFormatting = new List<OutputConfig> 
		{
		    new OutputConfig { outputType = OutputType.Control, index = 0, title="m", 
									scaleFactor=1 },
		    new OutputConfig { outputType = OutputType.State, index = 0, title="x0" },
		    new OutputConfig { outputType = OutputType.State, index = 1, title="x1" },
		};

	int dimensionX = 2;
	int dimensionM = 1;

	double dt = 0.1;
	int N = 51;

	var r = new Vector(dimensionX);
	var u = new Vector(dimensionM);
	var Q = new Matrix(dimensionX);
	var Z = new Matrix(dimensionM);

	r[0] = 0;
	r[1] = 0;

	u[0] = 0;

	Q[0, 0] = 1;
	Q[1, 1] = 1;

	Z[0, 0] = 1;

	var xInit = new Vector(dimensionX);
	xInit[0] = 1;

	double mu = 1;

	var functions = new CalcDelegate[]
		{
		    (k, deltaT, m, x) => x[1],
		    (k, deltaT, m, x) => mu * x[1] * (1 - x[0] * x[0]) - x[0] + m[0]
		};

	// Exact formulas for gradients calculation
	var gradientsA = new CalcDelegate[dimensionX, dimensionX];
	gradientsA[0, 0] = (k, deltaT, m, x) => 0;
	gradientsA[0, 1] = (k, deltaT, m, x) => 1;
	gradientsA[1, 0] = (k, deltaT, m, x) => -1 - 2 * x[0] * x[1];
	gradientsA[1, 1] = (k, deltaT, m, x) => 1 - x[0] * x[0];

	var gradientsB = new CalcDelegate[dimensionX, dimensionM];
	gradientsB[0, 0] = (k, deltaT, m, x) => 0;
	gradientsB[1, 0] = (k, deltaT, m, x) => 1;

	var dynamicSystem1 = SQ.Create(functions, gradientsA, gradientsB, Q, Z, r, u, xInit, dt, N,
				       (currCost, prevCost, iteration) => iteration > 20);
	dynamicSystem1.RunOptimization();
	
	dynamicSystem1.Output2Csv("VanDerPol_ExactGrads", lstFormatting);
	dynamicSystem1.OutputExecutionDetails(isExactGradient: true);

	// Numeric gradients calculation
	double delta = 0.1;
	var dynamicSystem2 = SQ.Create(functions, new Vector(dimensionM, delta), 
						  new Vector(dimensionX, delta),
					          Q, Z, r, u, xInit, dt, N,
					          (currCost, prevCost, iteration) => iteration > 20);
	dynamicSystem2.RunOptimization();

	dynamicSystem2.Output2Csv("VanDerPol_NumGrads", lstFormatting);
	dynamicSystem2.OutputExecutionDetails(isExactGradient: false);

	Console.WriteLine($"End \"{TITLE}\"{Environment.NewLine}{Environment.NewLine}");
}