车辆的动力学模型综合考虑了车辆运动过程中的惯性力和摩擦力的影响,其建模过程相对于运动学模型也稍微复杂一些。
1. 车体动力学建模
如下图所示,是车辆实际运动过程中的轮胎受力方向及运动的实际速度方向。
根据牛顿定律
(1)
(2)
横向加速度
(3)
摩擦力
(4)
(5)
为前轮绕重心自转产生的线速度
(6)
为后轮绕重心自转产生的线速度
(7)
由于摩擦角很小,采用近似计算
(8)
(9)
将(3),(4),(5),(8),(9)代入到(1)得到:
(10)
将(4),(5),(8),(9)代入到(2)得到:
(11)
故合并(10),(11)得到动力方程为:
(12)
2. 车体误差模型
车体误差模型解决的是,在速度恒定的情况,如何控制转向角 ,使车体按期望的角速度旋转,从而保证车体横向的控制误差最小。假设期望的角速度为 ,则期望的横向加速度为 。
(13)
(14)
横向位移加速度误差为:
(15)
车体的旋转角度误差为:
(16)
横向位移速度误差:
(17)
航向角速度误差:
(18)
航向角加速度误差:
(19)
车体旋转角加速度误差:
(20)
(21)
最终的误差状态方程为:
(22)
3. 横向LQR控制
令 , ,
则有:
在控制的某瞬时时刻,认为LQR控制器控制跟随的是轨迹的切线方向,也就是曲率0,对应的 和 为0。而期望轨迹的曲率变化可以转换为车辆控制的期望转向角 ,期望转向角 最终可以以前馈的方式加入到转向角输出控制中。
因此控制的状态方程为:
对状态方程进行离散化:
令:
,
则有:
至此,离散的状态方程建立完毕,接着和前面文章描述的LQR控制一样,套用公式即可。
迭代求解代数黎卡提方程:
当 ,其中eps为自己设定的精度值。如果达到该条件(设定的迭代精度和迭代次数)我们认为该方程已经收敛,此时根据获得的P矩阵求反馈控制增益K。
求取反馈增益矩
反馈控制量:
期望轨迹下,车辆期望的转角控制量为:
车辆最终的输出控制量为:
下面以Apollo的LQR横向控制代码为参考,分析LQR横向控制的代码实现:
误差模型的状态矩阵 的常量元素部分初始化。
Status LatController::Init(std::shared_ptr<DependencyInjector> injector,
const ControlConf *control_conf) {
control_conf_ = control_conf;
injector_ = injector;
if (!LoadControlConf(control_conf_)) {
AERROR << "failed to load control conf";
return Status(ErrorCode::CONTROL_COMPUTE_ERROR,
"failed to load control_conf");
}
// Matrix init operations.
const int matrix_size = basic_state_size_ + preview_window_;
matrix_a_ = Matrix::Zero(basic_state_size_, basic_state_size_);
matrix_ad_ = Matrix::Zero(basic_state_size_, basic_state_size_);
matrix_adc_ = Matrix::Zero(matrix_size, matrix_size);
/*
A matrix (Gear Drive)
[0.0, 1.0, 0.0, 0.0;
0.0, (-(c_f + c_r) / m) / v, (c_f + c_r) / m,
(l_r * c_r - l_f * c_f) / m / v;
0.0, 0.0, 0.0, 1.0;
0.0, ((lr * cr - lf * cf) / i_z) / v, (l_f * c_f - l_r * c_r) / i_z,
(-1.0 * (l_f^2 * c_f + l_r^2 * c_r) / i_z) / v;]
*/
matrix_a_(0, 1) = 1.0;
matrix_a_(1, 2) = (cf_ + cr_) / mass_;
matrix_a_(2, 3) = 1.0;
matrix_a_(3, 2) = (lf_ * cf_ - lr_ * cr_) / iz_;
matrix_a_coeff_ = Matrix::Zero(matrix_size, matrix_size);
matrix_a_coeff_(1, 1) = -(cf_ + cr_) / mass_;
matrix_a_coeff_(1, 3) = (lr_ * cr_ - lf_ * cf_) / mass_;
matrix_a_coeff_(3, 1) = (lr_ * cr_ - lf_ * cf_) / iz_;
matrix_a_coeff_(3, 3) = -1.0 * (lf_ * lf_ * cf_ + lr_ * lr_ * cr_) / iz_;
/*
b = [0.0, c_f / m, 0.0, l_f * c_f / i_z]^T
*/
matrix_b_ = Matrix::Zero(basic_state_size_, 1);
matrix_bd_ = Matrix::Zero(basic_state_size_, 1);
matrix_bdc_ = Matrix::Zero(matrix_size, 1);
matrix_b_(1, 0) = cf_ / mass_;
matrix_b_(3, 0) = lf_ * cf_ / iz_;
matrix_bd_ = matrix_b_ * ts_;
matrix_state_ = Matrix::Zero(matrix_size, 1);
matrix_k_ = Matrix::Zero(1, matrix_size);
matrix_r_ = Matrix::Identity(1, 1);
matrix_q_ = Matrix::Zero(matrix_size, matrix_size);
int q_param_size = control_conf_->lat_controller_conf().matrix_q_size();
int reverse_q_param_size =
control_conf_->lat_controller_conf().reverse_matrix_q_size();
if (matrix_size != q_param_size || matrix_size != reverse_q_param_size) {
const auto error_msg = absl::StrCat(
"lateral controller error: matrix_q size: ", q_param_size,
"lateral controller error: reverse_matrix_q size: ",
reverse_q_param_size,
" in parameter file not equal to matrix_size: ", matrix_size);
AERROR << error_msg;
return Status(ErrorCode::CONTROL_COMPUTE_ERROR, error_msg);
}
for (int i = 0; i < q_param_size; ++i) {
matrix_q_(i, i) = control_conf_->lat_controller_conf().matrix_q(i);
}
matrix_q_updated_ = matrix_q_;
InitializeFilters(control_conf_);
auto &lat_controller_conf = control_conf_->lat_controller_conf();
LoadLatGainScheduler(lat_controller_conf);
LogInitParameters();
enable_leadlag_ = control_conf_->lat_controller_conf()
.enable_reverse_leadlag_compensation();
if (enable_leadlag_) {
leadlag_controller_.Init(lat_controller_conf.reverse_leadlag_conf(), ts_);
}
enable_mrac_ =
control_conf_->lat_controller_conf().enable_steer_mrac_control();
if (enable_mrac_) {
mrac_controller_.Init(lat_controller_conf.steer_mrac_conf(),
vehicle_param_.steering_latency_param(), ts_);
}
enable_look_ahead_back_control_ =
control_conf_->lat_controller_conf().enable_look_ahead_back_control();
return Status::OK();
}
LQR横向控制主体计算函数:
Status LatController::ComputeControlCommand(
const localization::LocalizationEstimate *localization,
const canbus::Chassis *chassis,
const planning::ADCTrajectory *planning_published_trajectory,
ControlCommand *cmd) {
auto vehicle_state = injector_->vehicle_state();
auto target_tracking_trajectory = *planning_published_trajectory;
if (FLAGS_use_navigation_mode &&
FLAGS_enable_navigation_mode_position_update) {
auto time_stamp_diff =
planning_published_trajectory->header().timestamp_sec() -
current_trajectory_timestamp_;
auto curr_vehicle_x = localization->pose().position().x();
auto curr_vehicle_y = localization->pose().position().y();
double curr_vehicle_heading = 0.0;
const auto &orientation = localization->pose().orientation();
if (localization->pose().has_heading()) {
curr_vehicle_heading = localization->pose().heading();
} else {
curr_vehicle_heading =
common::math::QuaternionToHeading(orientation.qw(), orientation.qx(),
orientation.qy(), orientation.qz());
}
// new planning trajectory
if (time_stamp_diff > 1.0e-6) {
init_vehicle_x_ = curr_vehicle_x;
init_vehicle_y_ = curr_vehicle_y;
init_vehicle_heading_ = curr_vehicle_heading;
current_trajectory_timestamp_ =
planning_published_trajectory->header().timestamp_sec();
} else {
auto x_diff_map = curr_vehicle_x - init_vehicle_x_;
auto y_diff_map = curr_vehicle_y - init_vehicle_y_;
auto theta_diff = curr_vehicle_heading - init_vehicle_heading_;
auto cos_map_veh = std::cos(init_vehicle_heading_);
auto sin_map_veh = std::sin(init_vehicle_heading_);
auto x_diff_veh = cos_map_veh * x_diff_map + sin_map_veh * y_diff_map;
auto y_diff_veh = -sin_map_veh * x_diff_map + cos_map_veh * y_diff_map;
auto cos_theta_diff = std::cos(-theta_diff);
auto sin_theta_diff = std::sin(-theta_diff);
auto tx = -(cos_theta_diff * x_diff_veh - sin_theta_diff * y_diff_veh);
auto ty = -(sin_theta_diff * x_diff_veh + cos_theta_diff * y_diff_veh);
auto ptr_trajectory_points =
target_tracking_trajectory.mutable_trajectory_point();
std::for_each(
ptr_trajectory_points->begin(), ptr_trajectory_points->end(),
[&cos_theta_diff, &sin_theta_diff, &tx, &ty,
&theta_diff](common::TrajectoryPoint &p) {
auto x = p.path_point().x();
auto y = p.path_point().y();
auto theta = p.path_point().theta();
auto x_new = cos_theta_diff * x - sin_theta_diff * y + tx;
auto y_new = sin_theta_diff * x + cos_theta_diff * y + ty;
auto theta_new = common::math::NormalizeAngle(theta - theta_diff);
p.mutable_path_point()->set_x(x_new);
p.mutable_path_point()->set_y(y_new);
p.mutable_path_point()->set_theta(theta_new);
});
}
}
trajectory_analyzer_ =
std::move(TrajectoryAnalyzer(&target_tracking_trajectory));
// Transform the coordinate of the planning trajectory from the center of the
// rear-axis to the center of mass, if conditions matched
if (((FLAGS_trajectory_transform_to_com_reverse &&
vehicle_state->gear() == canbus::Chassis::GEAR_REVERSE) ||
(FLAGS_trajectory_transform_to_com_drive &&
vehicle_state->gear() == canbus::Chassis::GEAR_DRIVE)) &&
enable_look_ahead_back_control_) {
trajectory_analyzer_.TrajectoryTransformToCOM(lr_);
}
// Re-build the vehicle dynamic models at reverse driving (in particular,
// replace the lateral translational motion dynamics with the corresponding
// kinematic models)
if (vehicle_state->gear() == canbus::Chassis::GEAR_REVERSE) {
/*
A matrix (Gear Reverse)
[0.0, 0.0, 1.0 * v 0.0;
0.0, (-(c_f + c_r) / m) / v, (c_f + c_r) / m,
(l_r * c_r - l_f * c_f) / m / v;
0.0, 0.0, 0.0, 1.0;
0.0, ((lr * cr - lf * cf) / i_z) / v, (l_f * c_f - l_r * c_r) / i_z,
(-1.0 * (l_f^2 * c_f + l_r^2 * c_r) / i_z) / v;]
*/
cf_ = -control_conf_->lat_controller_conf().cf();
cr_ = -control_conf_->lat_controller_conf().cr();
matrix_a_(0, 1) = 0.0;
matrix_a_coeff_(0, 2) = 1.0;
} else {
/*
A matrix (Gear Drive)
[0.0, 1.0, 0.0, 0.0;
0.0, (-(c_f + c_r) / m) / v, (c_f + c_r) / m,
(l_r * c_r - l_f * c_f) / m / v;
0.0, 0.0, 0.0, 1.0;
0.0, ((lr * cr - lf * cf) / i_z) / v, (l_f * c_f - l_r * c_r) / i_z,
(-1.0 * (l_f^2 * c_f + l_r^2 * c_r) / i_z) / v;]
*/
cf_ = control_conf_->lat_controller_conf().cf();
cr_ = control_conf_->lat_controller_conf().cr();
matrix_a_(0, 1) = 1.0;
matrix_a_coeff_(0, 2) = 0.0;
}
matrix_a_(1, 2) = (cf_ + cr_) / mass_;
matrix_a_(3, 2) = (lf_ * cf_ - lr_ * cr_) / iz_;
matrix_a_coeff_(1, 1) = -(cf_ + cr_) / mass_;
matrix_a_coeff_(1, 3) = (lr_ * cr_ - lf_ * cf_) / mass_;
matrix_a_coeff_(3, 1) = (lr_ * cr_ - lf_ * cf_) / iz_;
matrix_a_coeff_(3, 3) = -1.0 * (lf_ * lf_ * cf_ + lr_ * lr_ * cr_) / iz_;
/*
b = [0.0, c_f / m, 0.0, l_f * c_f / i_z]^T
*/
matrix_b_(1, 0) = cf_ / mass_;
matrix_b_(3, 0) = lf_ * cf_ / iz_;
matrix_bd_ = matrix_b_ * ts_;
UpdateDrivingOrientation();
SimpleLateralDebug *debug = cmd->mutable_debug()->mutable_simple_lat_debug();
debug->Clear();
// Update state = [Lateral Error, Lateral Error Rate, Heading Error, Heading
// Error Rate, preview lateral error1 , preview lateral error2, ...]
UpdateState(debug);
UpdateMatrix();
// Compound discrete matrix with road preview model
UpdateMatrixCompound();
// Adjust matrix_q_updated when in reverse gear
int q_param_size = control_conf_->lat_controller_conf().matrix_q_size();
int reverse_q_param_size =
control_conf_->lat_controller_conf().reverse_matrix_q_size();
if (injector_->vehicle_state()->gear() == canbus::Chassis::GEAR_REVERSE) {
for (int i = 0; i < reverse_q_param_size; ++i) {
matrix_q_(i, i) =
control_conf_->lat_controller_conf().reverse_matrix_q(i);
}
} else {
for (int i = 0; i < q_param_size; ++i) {
matrix_q_(i, i) = control_conf_->lat_controller_conf().matrix_q(i);
}
}
// Add gain scheduler for higher speed steering
if (FLAGS_enable_gain_scheduler) {
matrix_q_updated_(0, 0) =
matrix_q_(0, 0) * lat_err_interpolation_->Interpolate(
std::fabs(vehicle_state->linear_velocity()));
matrix_q_updated_(2, 2) =
matrix_q_(2, 2) * heading_err_interpolation_->Interpolate(
std::fabs(vehicle_state->linear_velocity()));
common::math::SolveLQRProblem(matrix_adc_, matrix_bdc_, matrix_q_updated_,
matrix_r_, lqr_eps_, lqr_max_iteration_,
&matrix_k_);
} else {
common::math::SolveLQRProblem(matrix_adc_, matrix_bdc_, matrix_q_,
matrix_r_, lqr_eps_, lqr_max_iteration_,
&matrix_k_);
}
// feedback = - K * state
// Convert vehicle steer angle from rad to degree and then to steer degree
// then to 100% ratio
const double steer_angle_feedback = -(matrix_k_ * matrix_state_)(0, 0) * 180 /
M_PI * steer_ratio_ /
steer_single_direction_max_degree_ * 100;
//u = k*x
const double steer_angle_feedforward = ComputeFeedForward(debug->curvature());
//将期望轨迹曲率转换为车体期望的转向角
double steer_angle = 0.0;
double steer_angle_feedback_augment = 0.0;
// Augment the feedback control on lateral error at the desired speed domain
if (enable_leadlag_) {
if (FLAGS_enable_feedback_augment_on_high_speed ||
std::fabs(vehicle_state->linear_velocity()) < low_speed_bound_) {
steer_angle_feedback_augment =
leadlag_controller_.Control(-matrix_state_(0, 0), ts_) * 180 / M_PI *
steer_ratio_ / steer_single_direction_max_degree_ * 100;
if (std::fabs(vehicle_state->linear_velocity()) >
low_speed_bound_ - low_speed_window_) {
// Within the low-high speed transition window, linerly interplolate the
// augment control gain for "soft" control switch
steer_angle_feedback_augment = common::math::lerp(
steer_angle_feedback_augment, low_speed_bound_ - low_speed_window_,
0.0, low_speed_bound_, std::fabs(vehicle_state->linear_velocity()));
}
}
}
steer_angle = steer_angle_feedback + steer_angle_feedforward + //steer = ff + fb
steer_angle_feedback_augment;
// Compute the steering command limit with the given maximum lateral
// acceleration
const double steer_limit =
FLAGS_set_steer_limit ? std::atan(max_lat_acc_ * wheelbase_ /
(vehicle_state->linear_velocity() *
vehicle_state->linear_velocity())) *
steer_ratio_ * 180 / M_PI /
steer_single_direction_max_degree_ * 100
: 100.0;
const double steer_diff_with_max_rate =
FLAGS_enable_maximum_steer_rate_limit
? vehicle_param_.max_steer_angle_rate() * ts_ * 180 / M_PI /
steer_single_direction_max_degree_ * 100
: 100.0;
const double steering_position = chassis->steering_percentage();
// Re-compute the steering command if the MRAC control is enabled, with steer
// angle limitation and steer rate limitation
if (enable_mrac_) {
const int mrac_model_order = control_conf_->lat_controller_conf()
.steer_mrac_conf()
.mrac_model_order();
Matrix steer_state = Matrix::Zero(mrac_model_order, 1);
steer_state(0, 0) = chassis->steering_percentage();
if (mrac_model_order > 1) {
steer_state(1, 0) = (steering_position - pre_steering_position_) / ts_;
}
if (std::fabs(vehicle_state->linear_velocity()) >
control_conf_->minimum_speed_resolution()) {
mrac_controller_.SetStateAdaptionRate(1.0);
mrac_controller_.SetInputAdaptionRate(1.0);
} else {
mrac_controller_.SetStateAdaptionRate(0.0);
mrac_controller_.SetInputAdaptionRate(0.0);
}
steer_angle = mrac_controller_.Control(
steer_angle, steer_state, steer_limit, steer_diff_with_max_rate / ts_);
// Set the steer mrac debug message
MracDebug *mracdebug = debug->mutable_steer_mrac_debug();
Matrix steer_reference = mrac_controller_.CurrentReferenceState();
mracdebug->set_mrac_model_order(mrac_model_order);
for (int i = 0; i < mrac_model_order; ++i) {
mracdebug->add_mrac_reference_state(steer_reference(i, 0));
mracdebug->add_mrac_state_error(steer_state(i, 0) -
steer_reference(i, 0));
mracdebug->mutable_mrac_adaptive_gain()->add_state_adaptive_gain(
mrac_controller_.CurrentStateAdaptionGain()(i, 0));
}
mracdebug->mutable_mrac_adaptive_gain()->add_input_adaptive_gain(
mrac_controller_.CurrentInputAdaptionGain()(0, 0));
mracdebug->set_mrac_reference_saturation_status(
mrac_controller_.ReferenceSaturationStatus());
mracdebug->set_mrac_control_saturation_status(
mrac_controller_.ControlSaturationStatus());
}
pre_steering_position_ = steering_position;
debug->set_steer_mrac_enable_status(enable_mrac_);
// Clamp the steer angle with steer limitations at current speed
double steer_angle_limited =
common::math::Clamp(steer_angle, -steer_limit, steer_limit);
steer_angle = steer_angle_limited;
debug->set_steer_angle_limited(steer_angle_limited);
// Limit the steering command with the designed digital filter
steer_angle = digital_filter_.Filter(steer_angle);
steer_angle = common::math::Clamp(steer_angle, -100.0, 100.0);
// Check if the steer is locked and hence the previous steer angle should be
// executed
if (std::abs(vehicle_state->linear_velocity()) < FLAGS_lock_steer_speed &&
(vehicle_state->gear() == canbus::Chassis::GEAR_DRIVE ||
vehicle_state->gear() == canbus::Chassis::GEAR_REVERSE) &&
chassis->driving_mode() == canbus::Chassis::COMPLETE_AUTO_DRIVE) {
steer_angle = pre_steer_angle_;
}
// Set the steer commands
cmd->set_steering_target(common::math::Clamp(
steer_angle, pre_steer_angle_ - steer_diff_with_max_rate,
pre_steer_angle_ + steer_diff_with_max_rate));
cmd->set_steering_rate(FLAGS_steer_angle_rate);
pre_steer_angle_ = cmd->steering_target();
// compute extra information for logging and debugging
const double steer_angle_lateral_contribution =
-matrix_k_(0, 0) * matrix_state_(0, 0) * 180 / M_PI * steer_ratio_ /
steer_single_direction_max_degree_ * 100;
const double steer_angle_lateral_rate_contribution =
-matrix_k_(0, 1) * matrix_state_(1, 0) * 180 / M_PI * steer_ratio_ /
steer_single_direction_max_degree_ * 100;
const double steer_angle_heading_contribution =
-matrix_k_(0, 2) * matrix_state_(2, 0) * 180 / M_PI * steer_ratio_ /
steer_single_direction_max_degree_ * 100;
const double steer_angle_heading_rate_contribution =
-matrix_k_(0, 3) * matrix_state_(3, 0) * 180 / M_PI * steer_ratio_ /
steer_single_direction_max_degree_ * 100;
debug->set_heading(driving_orientation_);
debug->set_steer_angle(steer_angle);
debug->set_steer_angle_feedforward(steer_angle_feedforward);
debug->set_steer_angle_lateral_contribution(steer_angle_lateral_contribution);
debug->set_steer_angle_lateral_rate_contribution(
steer_angle_lateral_rate_contribution);
debug->set_steer_angle_heading_contribution(steer_angle_heading_contribution);
debug->set_steer_angle_heading_rate_contribution(
steer_angle_heading_rate_contribution);
debug->set_steer_angle_feedback(steer_angle_feedback);
debug->set_steer_angle_feedback_augment(steer_angle_feedback_augment);
debug->set_steering_position(steering_position);
debug->set_ref_speed(vehicle_state->linear_velocity());
ProcessLogs(debug, chassis);
return Status::OK();
}
更新车体当前的跟踪状态误差X:
void LatController::UpdateState(SimpleLateralDebug *debug) {
auto vehicle_state = injector_->vehicle_state();
if (FLAGS_use_navigation_mode) {
ComputeLateralErrors(
0.0, 0.0, driving_orientation_, vehicle_state->linear_velocity(),
vehicle_state->angular_velocity(), vehicle_state->linear_acceleration(),
trajectory_analyzer_, debug);
} else {
// Transform the coordinate of the vehicle states from the center of the
// rear-axis to the center of mass, if conditions matched
const auto &com = vehicle_state->ComputeCOMPosition(lr_);
ComputeLateralErrors(
com.x(), com.y(), driving_orientation_,
vehicle_state->linear_velocity(), vehicle_state->angular_velocity(),
vehicle_state->linear_acceleration(), trajectory_analyzer_, debug);
}
// State matrix update;
// First four elements are fixed;
if (enable_look_ahead_back_control_) {
matrix_state_(0, 0) = debug->lateral_error_feedback();
matrix_state_(2, 0) = debug->heading_error_feedback();
} else {
matrix_state_(0, 0) = debug->lateral_error();
matrix_state_(2, 0) = debug->heading_error();
}
matrix_state_(1, 0) = debug->lateral_error_rate();
matrix_state_(3, 0) = debug->heading_error_rate();
// Next elements are depending on preview window size;
for (int i = 0; i < preview_window_; ++i) {
const double preview_time = ts_ * (i + 1);
const auto preview_point =
trajectory_analyzer_.QueryNearestPointByRelativeTime(preview_time);
const auto matched_point = trajectory_analyzer_.QueryNearestPointByPosition(
preview_point.path_point().x(), preview_point.path_point().y());
const double dx =
preview_point.path_point().x() - matched_point.path_point().x();
const double dy =
preview_point.path_point().y() - matched_point.path_point().y();
const double cos_matched_theta =
std::cos(matched_point.path_point().theta());
const double sin_matched_theta =
std::sin(matched_point.path_point().theta());
const double preview_d_error =
cos_matched_theta * dy - sin_matched_theta * dx;
matrix_state_(basic_state_size_ + i, 0) = preview_d_error;
}
}
更新误差模型的状态矩阵A
void LatController::UpdateMatrix() {
double v;
// At reverse driving, replace the lateral translational motion dynamics with
// the corresponding kinematic models
if (injector_->vehicle_state()->gear() == canbus::Chassis::GEAR_REVERSE) {
v = std::min(injector_->vehicle_state()->linear_velocity(),
-minimum_speed_protection_);
matrix_a_(0, 2) = matrix_a_coeff_(0, 2) * v;
} else {
v = std::max(injector_->vehicle_state()->linear_velocity(),
minimum_speed_protection_);
matrix_a_(0, 2) = 0.0;
}
matrix_a_(1, 1) = matrix_a_coeff_(1, 1) / v;
matrix_a_(1, 3) = matrix_a_coeff_(1, 3) / v;
matrix_a_(3, 1) = matrix_a_coeff_(3, 1) / v;
matrix_a_(3, 3) = matrix_a_coeff_(3, 3) / v;
Matrix matrix_i = Matrix::Identity(matrix_a_.cols(), matrix_a_.cols());
matrix_ad_ = (matrix_i - ts_ * 0.5 * matrix_a_).inverse() *
(matrix_i + ts_ * 0.5 * matrix_a_); //离散状态矩阵A
}
LQR横向控制反馈增益矩阵K的计算函数,迭代代数黎卡提方程P,求取反馈增益矩阵K
void SolveLQRProblem(const Matrix &A, const Matrix &B, const Matrix &Q,
const Matrix &R, const Matrix &M, const double tolerance,
const uint max_num_iteration, Matrix *ptr_K) {
if (A.rows() != A.cols() || B.rows() != A.rows() || Q.rows() != Q.cols() ||
Q.rows() != A.rows() || R.rows() != R.cols() || R.rows() != B.cols() ||
M.rows() != Q.rows() || M.cols() != R.cols()) {
AERROR << "LQR solver: one or more matrices have incompatible dimensions.";
return;
}
Matrix AT = A.transpose();
Matrix BT = B.transpose();
Matrix MT = M.transpose();
// Solves a discrete-time Algebraic Riccati equation (DARE)
// Calculate Matrix Difference Riccati Equation, initialize P and Q
Matrix P = Q;
uint num_iteration = 0;
double diff = std::numeric_limits<double>::max();
while (num_iteration++ < max_num_iteration && diff > tolerance) {
Matrix P_next =
AT * P * A -
(AT * P * B + M) * (R + BT * P * B).inverse() * (BT * P * A + MT) + Q;
// check the difference between P and P_next
diff = fabs((P_next - P).maxCoeff());
P = P_next;
}
if (num_iteration >= max_num_iteration) {
ADEBUG << "LQR solver cannot converge to a solution, "
"last consecutive result diff is: "
<< diff;
} else {
ADEBUG << "LQR solver converged at iteration: " << num_iteration
<< ", max consecutive result diff.: " << diff;
}
*ptr_K = (R + BT * P * B).inverse() * (BT * P * A + MT);
}
评论(0)
您还未登录,请登录后发表或查看评论