面向阿克曼移动机器人(自行车模型)的LQR(最优二次型调节器)路径跟踪方法

        线性二次调节器(Linear Quadratic Regulator,LQR)是针对线性系统的最优控制方法。LQR 方法标准的求解体系是在考虑到损耗尽可能小的情况下, 以尽量小的代价平衡其他状态分量。一般情况下,线性系统在LQR 控制方法中用状态空间方程描述,性能能指标函数由二次型函数描述。

LQR 方法存在以下优点:

  1. 最小能量消耗和最高路径跟踪精度。
  2. 求解时能够考虑多状态情况。
  3. 鲁棒性较强。

缺点:

  1. 控制效果和模型精确程度有很大相关性。
  2. 实时计算状态反馈矩阵和控制增益。

一、系统模型

1.1 车辆模型

        一般来说阿克曼移动机器人可以简化为自行车模型,是一个非线性时变系统,工程上一般通过在平衡点附近差分线性化转化为线性系统来分析和控制,具体就不推导了,我直接给出模型。

1.2 线性系统状态反馈控制示意图

状态反馈是线性能控线性系统镇定的一个有效方法,主要是通过极点配置方法寻找一组非正的闭环极点使得闭环系统大范围渐进稳定。

 A,B,C分别代表系统矩阵、输入矩阵和输出矩阵,K是待设计的状态反馈增益。

二、控制器设计

2.1 代价函数泛函设计

最优控制里,代价函数一般设计为控制性能和控制代价的范数加权和,形式如下

其中,期望和实际的误差系统定义为

2.2 最优状态反馈控制律推导

当想要状态与期望状态之间的误差越差越小,同时控制消耗更少的能量。求解极小值点时,新定义拉格朗日函数如下 

在拉格朗日函数基础上对各个优化变量的一阶导为零 ,得

 当\Delta_{x(t)}=0时候,

推导LQR控制律时候,设\lambda(t)=P_tx(t) ,求偏导后可得

 得

 由于状态量初始不为零,只能是

又由于当上述方程成立时候,P收敛到了期望的范围 ,\dot{P}为零,因此得到立卡提方程形式的矩阵微分方程

综上,通过迭代或者近似方法求解上述立卡提方程后,最优的控制律为

 2.3 连续立卡提方程求解流程

For~~ iteration=1: iteration_{max}\\ ~~~~~~~\bar{P}=Q+A^{\top}PB(R+B^TPB)^{-1}B^{\top}PA\\ ~~~~~~~criteria=\overline\lambda(\vert\bar{P}-P\vert)\\ ~~~~~~~IF~~criteria<\epsilon\\ ~~~~~~~~~~Break\\ ~~~~~~P=\bar{P}

三、具体实现代码

3.1 main函数

close all
clear;
clc; 
cx = [];
cy= [];
y0 = @(t_step) 10*sin(2 * t_step + 1);
x0_dot= @(t_step) 5 * 2 * cos(2 * t_step + 1);
x0 = @(t_step) -40*cos(t_step + 0.5);
for theta=0:pi/200:2*pi
    cx(end + 1) = x0(theta);
    cy(end + 1) = y0(theta);
end
refer_path_primary= [cx', cy'];
x = refer_path_primary(:, 1)';
y = refer_path_primary(:, 2)';
points = [x; y]';
ds = 0.1 ;%等距插值处理的间隔
distance = [0, cumsum(hypot(diff(x, 1), diff(y, 1)))]';
distance_specific = 0:ds:distance(end);
hypot(diff(x, 1), diff(y, 1));
diff(x, 1);
diff(y, 1);
s = 0:ds:distance(end);
refer_path= interp1(distance, points, distance_specific, 'spline');
figure(1)
% 绘制拟合曲线
plot(refer_path(:, 1), refer_path(:, 2),'b-',x, y,'r.');
hold on
refer_path_x = refer_path(:,1);  % x
refer_path_y = refer_path(:,2); % y
for i=1:length(refer_path)
    if i==1
         dx = refer_path(i + 1, 1) - refer_path(i, 1);
         dy = refer_path(i + 1, 2) - refer_path(i, 2);
         ddx = refer_path(3, 1) + refer_path(1, 1) - 2 * refer_path(2, 1);
         ddy = refer_path(3, 2) + refer_path(1, 2) - 2 * refer_path(2, 2);
    elseif  i==length(refer_path)
         dx = refer_path(i, 1) - refer_path(i - 1, 1);
         dy = refer_path(i, 2) - refer_path(i - 1, 2);
         ddx = refer_path(i, 1) + refer_path(i - 2, 1) - 2 * refer_path(i - 1, 1);
         ddy = refer_path(i, 2) + refer_path(i - 2, 2) - 2 * refer_path(i - 1, 2);
    else
          dx = refer_path(i + 1, 1) - refer_path(i, 1);
          dy = refer_path(i + 1, 2) - refer_path(i, 2);
          ddx = refer_path(i + 1, 1) + refer_path(i - 1, 1) - 2 * refer_path(i, 1);
          ddy = refer_path(i + 1, 2) + refer_path(i - 1, 2) - 2 * refer_path(i, 2);
    end
    refer_path(i,3)=atan2(dy, dx);%yaw
    refer_path(i,4)=(ddy * dx - ddx * dy) / ((dx ^ 2 + dy ^ 2) ^ (3 / 2));
end
figure(2)
plot(refer_path(:, 3),'b-');
figure(3)
plot(refer_path(:, 4),'b-')
% 
%%目标及初始状态
L=2;%车辆轴距
v=2;%初始速度
dt=0.05;%时间间隔
goal=refer_path(end,1:2);
x_0=refer_path_x(1);
y_0=refer_path_y(1);
psi_0 = refer_path(1, 3);
% %运动学模型
ugv = KinematicModel(x_0, y_0, psi_0, v, dt, L);
Q = eye(3) * 3.0;
R = eye(2) * 2.0;
robot_state = zeros(4, 1);
step_points = length(refer_path(:, 1));
for i=1:1:step_points
    robot_state(1)=ugv.x;
    robot_state(2)=ugv.y;
    robot_state(3)=ugv.psi;
    robot_state(4)=ugv.v;
    [e, k, ref_yaw, min_idx] = calc_track_error(robot_state(1), robot_state(2), refer_path);
    ref_delta = atan2(L * k, 1);
    [A, B] = state_space( robot_state(4), ref_delta, ref_yaw, dt, L);
    delta = lqr(robot_state, refer_path, min_idx, A, B, Q, R);
    delta = delta + ref_delta;
    [ugv.x, ugv.y, ugv.psi, ugv.v] = update(0, delta, dt, L, robot_state(1), robot_state(2),robot_state(3), robot_state(4));
    ugv.record_x(end + 1 ) = ugv.x;
    ugv.record_y(end + 1 ) = ugv.y;
    ugv.record_psi(end + 1 ) = ugv.psi;
     ugv.record_phy(end + 1 ) = ref_delta;
end
figure(4)
% 绘制拟合曲线
%  scr_size = get(0,'screensize');
% set(gcf,'outerposition', [1 1 scr_size(4) scr_size(4)]);
plot(ugv.record_x , ugv.record_y, Color='m',LineStyle='--',LineWidth=2);
axis([-40,40,-40,40])
grid on
hold on
% 绘制车辆曲线
axis equal
for ii = 1:1:length(ugv.record_x)
       
        h =  PlotCarbox(ugv.record_x(ii), ugv.record_y(ii), ugv.record_psi(ii), 'Color', 'r',LineWidth=2);
        h1 = plot(ugv.record_x(1:ii), ugv.record_y(1:ii),'Color', 'b');
        th1 = text(ugv.record_x(ii), ugv.record_y(ii)+10, ['#car', num2str(1)], 'Color', 'm');
        set(th1, {'HorizontalAlignment'},{'center'});
        h2 = PlotCarWheels(ugv.record_x(ii), ugv.record_y(ii), ugv.record_psi(ii),ugv.record_phy(ii),'k',LineWidth=2);
        h3 = plot(ugv.record_x(1:ii) , ugv.record_y(1:ii), Color='b',LineStyle='-',LineWidth=4);
        drawnow
        delete(h); delete(h1);delete(th1);delete(h3);
        for jj = 1:1:size(h2)
            delete(h2{jj});
        end
end
% 
function [P] =  cal_Ricatti(A, B, Q, R)
 Qf = Q;
 P = Qf;
  iter_max = 100;
  Eps = 1e-4;
 for step = 1:1:iter_max
     P_bar = Q + A' * P * A - A' * P * B * pinv(R + B' * P *B) * B' * P * A;
     criteria = max(abs(P_bar - P));
     if  criteria < Eps
         break;
     end
     P = P_bar;
 end

end
%%LQR控制器
function[u_star]=lqr(robot_state, refer_path, s0, A, B, Q, R)
    x = robot_state(1:3) - refer_path(s0,1:3)';
    P =  cal_Ricatti(A, B, Q, R);
    K= -pinv(R + B' * P * B) * B' * P * A;
    u = K * x;%状态反馈
    u_star = u(2);
end

function [e, k, yaw, min_idx]=calc_track_error(x, y, refer_path)
    p_num = length(refer_path);
    d_x = zeros(p_num, 1);
    d_y = zeros(p_num, 1);
    d = zeros(p_num, 1);
    for i=1:1:p_num 
        d_x(i) = refer_path(i, 1) - x; 
        d_y(i) = refer_path(i, 2) - y;
    end
    for i=1:1:p_num 
        d(i) = sqrt(d_x(i) ^2 + d_y(i) ^ 2) ;
    end
    [~, min_idx] = min(d); 
    yaw = refer_path(min_idx, 3);
    k = refer_path(min_idx, 4);
    angle= normalize_angle(yaw - atan2(d_y(min_idx), d_x(min_idx)));
    e = d(min_idx);
    if angle < 0
        e = e*(-1);
    end
end
%%将角度取值范围限定为[-pi,pi]
function [angle]=normalize_angle(angle)
    while angle > pi
        angle = angle - 2*pi;
    end
    while angle < pi
        angle = angle + 2*pi;
    end
end
function [x_next, y_next, psi_next, v_next] = update(a, delta_f, dt, L, x, y, psi, v)
     x_next = x + v * cos(psi) * dt;
     y_next = y + v * sin(psi) * dt;
     psi_next = psi + v / L * tan(delta_f) * dt;
     v_next = v + a * dt;
end
function [A, B]=state_space(v, ref_delta, ref_yaw, dt, L)
    A=[ 1.0, 0.0,  -v * dt * sin(ref_yaw);
           0.0, 1.0,   v *  dt * cos(ref_yaw);
           0.0, 0.0,   1.0  ];
    B =[ dt * cos(ref_yaw), 0;
           dt * sin(ref_yaw), 0;
           dt * tan(ref_delta) / L, v * dt / (L * cos(ref_delta) * cos(ref_delta))];
end
function h = PlotCarbox(x, y, theta, varargin)
Params = GetVehicleParams();

carbox = [-Params.Lr -Params.Lb/2; Params.Lw+Params.Lf -Params.Lb/2; Params.Lw+Params.Lf Params.Lb/2; -Params.Lr Params.Lb/2];
carbox = [carbox; carbox(1, :)];

transformed_carbox = [carbox ones(5, 1)] * GetTransformMatrix(x, y, theta)';
h = plot(transformed_carbox(:, 1), transformed_carbox(:, 2), varargin{:});

end
function hs = PlotCarWheels(x, y, theta, phy, varargin)
Params = GetVehicleParams();

wheel_box = [-Params.wheel_radius -Params.wheel_width / 2;
    +Params.wheel_radius -Params.wheel_width / 2;
    +Params.wheel_radius +Params.wheel_width / 2;
    -Params.wheel_radius +Params.wheel_width / 2];

front_x = x + Params.Lw * cos(theta);
front_y = y + Params.Lw * sin(theta);
track_width_2 = (Params.Lb - Params.wheel_width / 2) / 2;

boxes = {
    TransformBox(wheel_box, x - track_width_2 * sin(theta), y + track_width_2 * cos(theta), theta);
    TransformBox(wheel_box, x + track_width_2 * sin(theta), y - track_width_2 * cos(theta), theta);
    TransformBox(wheel_box, front_x - track_width_2 * sin(theta), front_y + track_width_2 * cos(theta), theta+phy);
    TransformBox(wheel_box, front_x + track_width_2 * sin(theta), front_y - track_width_2 * cos(theta), theta+phy);
};

hs = cell(4, 1);
for ii = 1:4
    hs{ii} = fill(boxes{ii}(:, 1), boxes{ii}(:, 2), varargin{:});
end

end

function transformed = TransformBox(box, x, y, theta)
transformed = [box; box(1, :)];
transformed = [transformed ones(5, 1)] * GetTransformMatrix(x, y, theta)';
transformed = transformed(:, 1:2);
end
function mat = GetTransformMatrix(x, y, theta)
mat = [ ...
    cos(theta) -sin(theta), x; ...
    sin(theta) cos(theta), y; ...
    0 0 1];
end

3.2 运动学结构体: 

classdef  KinematicModel<handle
    properties
        x;
        y;
        psi;
        v;
        dt;
        L;
        record_x;
        record_y;
        record_psi;
        record_phy;
    end
    methods
        function self=KinematicModel(x, y, psi, v, dt, L)
            self.x=x;
            self.y=y;
            self.psi=psi;
            self.v = v;
            self.L = L;
            % 实现是离散的模型
            self.dt = dt;
            self.record_x = [];
            self.record_y= [];
            self.record_psi= [];
            self.record_phy= [];
        end
end
end

四、仿真参数和效果

 4.1 参数配置

%期望轨迹
y0 = @(t_step) 10*sin(2 * t_step + 1);
x0_dot= @(t_step) 5 * 2 * cos(2 * t_step + 1);
L=2;%车辆轴距
v=2;%初始速度
dt=0.05;%时间间隔
Q = eye(3) * 3.0;
R = eye(2) * 2.0;
robot_state = zeros(4, 1);
VehicleParams.Lw = 2.8 * 2; % wheel base
VehicleParams.Lf = 0.96 * 2; % front hang
VehicleParams.Lr = 0.929 * 2; % rear hang
VehicleParams.Lb = 1.942 * 2; % width
VehicleParams.Ll = VehicleParams.Lw + VehicleParams.Lf + VehicleParams.Lr; % length
VehicleParams.f2x = 1/4 * (3*VehicleParams.Lw + 3*VehicleParams.Lf - VehicleParams.Lr);
VehicleParams.r2x = 1/4 * (VehicleParams.Lw + VehicleParams.Lf - 3*VehicleParams.Lr);
VehicleParams.radius = 1/2 * sqrt((VehicleParams.Lw + VehicleParams.Lf + VehicleParams.Lr) ^ 2 / 4 + VehicleParams.Lb ^ 2);
VehicleParams.a_max = 0.5;
VehicleParams.v_max = 2.5;
VehicleParams.phi_max = 0.7;
VehicleParams.omega_max = 0.5;
% for wheel visualization
VehicleParams.wheel_radius = 0.32*2;
VehicleParams.wheel_width = 0.22*2;
iter_max = 100;
Eps = 1e-4;

  4.1 仿真效果

本文来自互联网用户投稿,该文观点仅代表作者本人,不代表本站立场。本站仅提供信息存储空间服务,不拥有所有权,不承担相关法律责任。如若转载,请注明出处:/a/762476.html

如若内容造成侵权/违法违规/事实不符,请联系我们进行投诉反馈qq邮箱809451989@qq.com,一经查实,立即删除!

相关文章

汇聚荣拼多多电商好不好?

拼多多电商好不好?这是一个值得探讨的问题。拼多多作为中国领先的电商平台之一&#xff0c;以其独特的商业模式和创新的营销策略吸引了大量用户。然而&#xff0c;对于这个问题的回答并不是简单的好或不好&#xff0c;而是需要从多个方面进行综合分析。 一、商品质量 来看拼多…

【源码+文档+调试讲解】居家养老系统

摘要 随着信息技术在管理上越来越深入而广泛的应用&#xff0c;管理信息系统的实施在技术上已逐步成熟。本文介绍了居家养老系统的开发全过程。通过分析高校学生综合素质评价管理方面的不足&#xff0c;创建了一个计算机管理居家养老系统的方案。文章介绍了居家养老系统的系统分…

jvm性能监控常用工具

在java的/bin目录下有许多java自带的工具。 我们常用的有 基础工具 jar:创建和管理jar文件 java&#xff1a;java运行工具&#xff0c;用于运行class文件或jar文件 javac&#xff1a;java的编译器 javadoc&#xff1a;java的API文档生成工具 性能监控和故障处理 jps jstat…

Sourcecodester Fantastic Blog CMS v1.0 SQL 注入漏洞(CVE-2022-28512)

前言 CVE-2022-28512 是一个存在于 Sourcecodester Fantastic Blog CMS v1.0 中的 SQL 注入漏洞。攻击者可以通过 "/fantasticblog/single.php" 中的 id 参数注入恶意 SQL 查询&#xff0c;从而获得对数据库的未经授权的访问和控制。 漏洞详细信息 漏洞描述: 该漏…

JavaScript将参数传递给事件处理程序

本篇文件我们将实现导航栏中&#xff0c;选中时候&#xff0c;会将您选中的进行高亮显示&#xff1b; ● 首先我们来获取我们想要的HTML元素 const nav document.querySelector(.nav);● 接着我们来写选中的高亮显示 nav.addEventListener(mouseover, function (e) { //鼠…

内网穿透小工具

内网穿透小工具 前言 当在本地或者虚拟机&#xff0c;内网搭建了项目&#xff0c;数据库。可是在外网无法访问。下面的两款小工具可以暂时实现内网穿透能力。&#xff08;不支持自定义域名&#xff0c;但是不限制隧道数量&#xff01;且免费&#xff01;免费&#xff01;免费…

【小贪】项目实战——Zero-shot根据文字提示分割出图片目标掩码

目标描述 给定RGB视频或图片&#xff0c;目标是分割出图像中的指定目标掩码。我们需要复现两个Zero-shot的开源项目&#xff0c;分别为IDEA研究院的GroundingDINO和Facebook的SAM。首先使用目标检测方法GroundingDINO&#xff0c;输入想检测目标的文字提示&#xff0c;可以获得…

互联网框架五层模型详解

注&#xff1a;机翻&#xff0c;未校对。 What is the Five Layers Model? The Framework of the Internet Explained 五层模型互联网框架解释 Computer Networks are a beautiful, amazing topic. Networks involve so much knowledge from different fields, from physics…

[OHOS_ERROR]: Please call hb utilities inside ohos source directory

当执行hb set报如下错误时&#xff1a;原因时重新拉取了源码&#xff0c;且源码路径被改了 [OHOS_ERROR]: Please call hb utilities inside ohos source directory 【解决办法】 卸载hb并在源码路径下重新安装 python3 -m pip uninstall ohos-build 安装hb python3 -m pi…

python-逻辑语句

if else语句 不同于C&#xff1a;else if range语句&#xff1a; continue continue的作用是&#xff1a; 中断所在循环的当次执行&#xff0c;直接进入下一次 continue在嵌套循环中的应用 break 直接结束所在的循环 break在嵌套循环中的应用 continue和break&#xff0c;在…

力扣:LCR 024. 反转链表(Java)

目录 题目描述&#xff1a;示例 1&#xff1a;示例 2&#xff1a;代码实现&#xff1a; 题目描述&#xff1a; 给定单链表的头节点 head &#xff0c;请反转链表&#xff0c;并返回反转后的链表的头节点。 示例 1&#xff1a; 输入&#xff1a;head [1,2,3,4,5] 输出&#x…

【嵌入式DIY实例】- LCD ST7735显示DHT11传感器数据

LCD ST7735显示DHT11传感器数据 文章目录 LCD ST7735显示DHT11传感器数据1、硬件准备与接线2、代码实现本文介绍如何将 ESP8266 NodeMCU 板 (ESP-12E) 与 DHT11 (RHT01) 数字湿度和温度传感器连接。 NodeMCU 从 DHT11 传感器读取温度(以 C 为单位)和湿度(以 rH% 为单位)值,…

1.5 Canal 数据同步工具详细教程

欢迎来到我的博客&#xff0c;很高兴能够在这里和您见面&#xff01;欢迎订阅相关专栏&#xff1a; ⭐️ 全网最全IT互联网公司面试宝典&#xff1a;收集整理全网各大IT互联网公司技术、项目、HR面试真题. ⭐️ AIGC时代的创新与未来&#xff1a;详细讲解AIGC的概念、核心技术、…

【你也能从零基础学会网站开发】关系型数据库中的表(Table)设计结构以及核心组成部分

&#x1f680; 个人主页 极客小俊 ✍&#x1f3fb; 作者简介&#xff1a;程序猿、设计师、技术分享 &#x1f40b; 希望大家多多支持, 我们一起学习和进步&#xff01; &#x1f3c5; 欢迎评论 ❤️点赞&#x1f4ac;评论 &#x1f4c2;收藏 &#x1f4c2;加关注 关系型数据库中…

FTP 文件传输协议:概念、工作原理;上传下载操作步骤

目录 FTP 概念 工作原理 匿名用户 授权用户 FTP软件包 匿名用户上传下载实验步骤 环境配置 下载 上传 wget 授权用户上传下载步骤 root用户登录FTP步骤 监听 设置端口号范围 修改用户家目录 匿名用户 授权用户 FTP 概念 FTP&#xff08;File Transfer Prot…

C语言之线程的学习

线程属于某一个进程 共同点&#xff1a;都能并发 线程共享变量&#xff0c;进程不共享。 多线程任务中&#xff0c;其中某一个线程调用了exit了&#xff0c;其他线程会跟着一起退出 如果是特定的线程就调用pthread_exit 失败返回的是错误号 下面也是

VSCode无法识别 node、npm

一、前提 电脑新安装了node.js&#xff0c;在cmd查看node和npm版本没有问题&#xff0c;但是在VSCode无法识别 1.cmd查看版本&#xff1a; 2.VSCode报错信息&#xff1a; 无法将“npm”项识别为 cmdlet、函数、脚本文件或可运行程序的名称。请检查名称的拼写&#xff0c;如果…

C#——Property属性详情

属性 属性&#xff08;Property&#xff09;是类&#xff08;class&#xff09;、结构体&#xff08;structure&#xff09;和接口&#xff08;interface&#xff09;的成员&#xff0c;类或结构体中的成员变量称为字段&#xff0c;属性是字段的扩展&#xff0c;使用访问器&am…

视频转音频:怎样提取视频中的音频?6个提取音频的小技巧(建议收藏)

怎样提取视频中的音频&#xff1f;当我们想从视频中提取出声音时&#xff0c;通常会遇到很多问题。无论是想单独提取出视频里的音频&#xff0c;还是把它转成方便储存或者分享的音频格式&#xff0c;这都会涉及到视频转音频的一个需求。因此&#xff0c;在这篇指南里&#xff0…

如何编写高质量更优雅的代码(Java)

1、函数式接口—FunctionalInterface 好处&#xff1a;高逼格、代码收拢、解藕、统一处理 适用范围&#xff1a;具有共性的接口调用代码 举个栗子&#xff1a; 在我们平时的微服务开发中&#xff0c;调用其他服务的接口&#xff0c;通常要把接口调用部分做异常处理(try catch…