I have a question on determining the variance in the Econometrics Toolbox example "Using Bootstrapping and Filtered Historical Simulation to Evaluate Market Risk". Everything is pretty straitforward till the end of the example where they use the filter function to determine the portfolio returns. Here is the command.
portfolioReturns = filter(fit, bootstrappedResiduals, …
‘Y0’, Y0, ‘Z0’, Z0, ‘V0’, V0);
I would like to obtain the conditional variance paths. Hence I used the command
[V, Y] = filter(fit, bootstrappedResiduals, …
‘Y0’, Y0, ‘Z0’, Z0, ‘V0’, V0);
Here is a plot of the first 10 horizon results for V

Here is the plot of the firt 10 horizon results for Y

The results for Y look like one would expect – plus and minus values around zero. However, the results for V should all be positive since according to the filter function for conditional variances V is the filtered variance result. When comparing the value of V to that of the variable portfolioReturns I fould them to be identical.

I can always square the value of the portfolioReturns to get the squarded returns. Is this the value of the conditional variance for the FHS simulation?I have a question on determining the variance in the Econometrics Toolbox example "Using Bootstrapping and Filtered Historical Simulation to Evaluate Market Risk". Everything is pretty straitforward till the end of the example where they use the filter function to determine the portfolio returns. Here is the command.
portfolioReturns = filter(fit, bootstrappedResiduals, …
‘Y0’, Y0, ‘Z0’, Z0, ‘V0’, V0);
I would like to obtain the conditional variance paths. Hence I used the command
[V, Y] = filter(fit, bootstrappedResiduals, …
‘Y0’, Y0, ‘Z0’, Z0, ‘V0’, V0);
Here is a plot of the first 10 horizon results for V

Here is the plot of the firt 10 horizon results for Y

The results for Y look like one would expect – plus and minus values around zero. However, the results for V should all be positive since according to the filter function for conditional variances V is the filtered variance result. When comparing the value of V to that of the variable portfolioReturns I fould them to be identical.

I can always square the value of the portfolioReturns to get the squarded returns. Is this the value of the conditional variance for the FHS simulation? I have a question on determining the variance in the Econometrics Toolbox example "Using Bootstrapping and Filtered Historical Simulation to Evaluate Market Risk". Everything is pretty straitforward till the end of the example where they use the filter function to determine the portfolio returns. Here is the command.
portfolioReturns = filter(fit, bootstrappedResiduals, …
‘Y0’, Y0, ‘Z0’, Z0, ‘V0’, V0);
I would like to obtain the conditional variance paths. Hence I used the command
[V, Y] = filter(fit, bootstrappedResiduals, …
‘Y0’, Y0, ‘Z0’, Z0, ‘V0’, V0);
Here is a plot of the first 10 horizon results for V

Here is the plot of the firt 10 horizon results for Y

The results for Y look like one would expect – plus and minus values around zero. However, the results for V should all be positive since according to the filter function for conditional variances V is the filtered variance result. When comparing the value of V to that of the variable portfolioReturns I fould them to be identical.

I can always square the value of the portfolioReturns to get the squarded returns. Is this the value of the conditional variance for the FHS simulation? conditional variance, filtered historical simulation MATLAB Answers — New Questions

​

I have used application compiler in Matlab2023b to pack two executable applications (e.g. App1.exe, App2.exe).
I also chose the setting "Do not display the Windows Command Shell (console) for execution" in my project file for both of them.
App2.exe will be called inside one of the function in App1.exe, which is hard coded.
I want to run App1.exe in a batch mode and i also wrote a batch file with one command to run it.
In the terminal the console messages of App1.exe don’t show up (which is expected), but the console messages of App2.exe still exist.
Actually i also tried in the batch file like "App1.exe >nul 2>&1" which makes the situation even worse, console messages of App1.exe also shou up after this adjustment. And i also already use the @echo off, but i think it was not the solution.
How can i also hide the console messages of the App2.exe?
btw. this problem didn’t show in Matlab202b before.I have used application compiler in Matlab2023b to pack two executable applications (e.g. App1.exe, App2.exe).
I also chose the setting "Do not display the Windows Command Shell (console) for execution" in my project file for both of them.
App2.exe will be called inside one of the function in App1.exe, which is hard coded.
I want to run App1.exe in a batch mode and i also wrote a batch file with one command to run it.
In the terminal the console messages of App1.exe don’t show up (which is expected), but the console messages of App2.exe still exist.
Actually i also tried in the batch file like "App1.exe >nul 2>&1" which makes the situation even worse, console messages of App1.exe also shou up after this adjustment. And i also already use the @echo off, but i think it was not the solution.
How can i also hide the console messages of the App2.exe?
btw. this problem didn’t show in Matlab202b before. I have used application compiler in Matlab2023b to pack two executable applications (e.g. App1.exe, App2.exe).
I also chose the setting "Do not display the Windows Command Shell (console) for execution" in my project file for both of them.
App2.exe will be called inside one of the function in App1.exe, which is hard coded.
I want to run App1.exe in a batch mode and i also wrote a batch file with one command to run it.
In the terminal the console messages of App1.exe don’t show up (which is expected), but the console messages of App2.exe still exist.
Actually i also tried in the batch file like "App1.exe >nul 2>&1" which makes the situation even worse, console messages of App1.exe also shou up after this adjustment. And i also already use the @echo off, but i think it was not the solution.
How can i also hide the console messages of the App2.exe?
btw. this problem didn’t show in Matlab202b before. batch, executable MATLAB Answers — New Questions

​

Hi, I have a uitable and have a copy function as below. This works fine if the table is solely numbers

function CopyTableData(app,src,event)
tbl = event.ContextObject;
% fieldnames(event)
% event.InteractionInformation
s=tbl.Data;
size_s = size(s);

str = ”;
for i=1:size_s(1)
for j=1:size_s(2)
str = sprintf ( ‘%s%ft’, str, s(i,j) );
end
str = sprintf ( ‘%sn’, str );
end
clipboard(‘copy’,str);
end

There are instances where my table will be a cell array and I have modified the function (for the str = sprintf ( ‘%s%ft’, str, s{i,j}) to:
function CopyTableData(app,src,event)
tbl = event.ContextObject;
% fieldnames(event)
% event.InteractionInformation
s=tbl.Data;
size_s = size(s);

str = ”;
for i=1:size_s(1)
for j=1:size_s(2)
try
str = sprintf ( ‘%s%ft’, str, s(i,j) );
catch
str = sprintf ( ‘%s%ft’, str, s{i,j} );
end
end
str = sprintf ( ‘%sn’, str );
end
clipboard(‘copy’,str);
end

However, when my data is as below (cellarray) in which the contents of all columns in the cellarray are numbers, but the very last column ‘info" are charss

I get this when pasted into excel (notice the last column is incorrect and actually split into 2)

this is how I created my dummy data:

uit=app.UITableSummary;
dx=9437;
z1=5795;
z2=5775.2
dz=z2-z1
tilt_mrads=-2.1084
A={dx,z1,z2,dz,tilt_mrads,’info’};
d=uit.Data; % Get Current table data
uit.Data = [d; A]; %Add new row to current data

I kept on adding this and just manually replaced the ‘info’ value in the last column (’10’, ’12’, etc..) and the numeric values in the 2nd from last column (-2.11, -1.59, etc..)

thanks for any helpHi, I have a uitable and have a copy function as below. This works fine if the table is solely numbers

function CopyTableData(app,src,event)
tbl = event.ContextObject;
% fieldnames(event)
% event.InteractionInformation
s=tbl.Data;
size_s = size(s);

str = ”;
for i=1:size_s(1)
for j=1:size_s(2)
str = sprintf ( ‘%s%ft’, str, s(i,j) );
end
str = sprintf ( ‘%sn’, str );
end
clipboard(‘copy’,str);
end

There are instances where my table will be a cell array and I have modified the function (for the str = sprintf ( ‘%s%ft’, str, s{i,j}) to:
function CopyTableData(app,src,event)
tbl = event.ContextObject;
% fieldnames(event)
% event.InteractionInformation
s=tbl.Data;
size_s = size(s);

str = ”;
for i=1:size_s(1)
for j=1:size_s(2)
try
str = sprintf ( ‘%s%ft’, str, s(i,j) );
catch
str = sprintf ( ‘%s%ft’, str, s{i,j} );
end
end
str = sprintf ( ‘%sn’, str );
end
clipboard(‘copy’,str);
end

However, when my data is as below (cellarray) in which the contents of all columns in the cellarray are numbers, but the very last column ‘info" are charss

I get this when pasted into excel (notice the last column is incorrect and actually split into 2)

this is how I created my dummy data:

uit=app.UITableSummary;
dx=9437;
z1=5795;
z2=5775.2
dz=z2-z1
tilt_mrads=-2.1084
A={dx,z1,z2,dz,tilt_mrads,’info’};
d=uit.Data; % Get Current table data
uit.Data = [d; A]; %Add new row to current data

I kept on adding this and just manually replaced the ‘info’ value in the last column (’10’, ’12’, etc..) and the numeric values in the 2nd from last column (-2.11, -1.59, etc..)

thanks for any help Hi, I have a uitable and have a copy function as below. This works fine if the table is solely numbers

function CopyTableData(app,src,event)
tbl = event.ContextObject;
% fieldnames(event)
% event.InteractionInformation
s=tbl.Data;
size_s = size(s);

str = ”;
for i=1:size_s(1)
for j=1:size_s(2)
str = sprintf ( ‘%s%ft’, str, s(i,j) );
end
str = sprintf ( ‘%sn’, str );
end
clipboard(‘copy’,str);
end

There are instances where my table will be a cell array and I have modified the function (for the str = sprintf ( ‘%s%ft’, str, s{i,j}) to:
function CopyTableData(app,src,event)
tbl = event.ContextObject;
% fieldnames(event)
% event.InteractionInformation
s=tbl.Data;
size_s = size(s);

str = ”;
for i=1:size_s(1)
for j=1:size_s(2)
try
str = sprintf ( ‘%s%ft’, str, s(i,j) );
catch
str = sprintf ( ‘%s%ft’, str, s{i,j} );
end
end
str = sprintf ( ‘%sn’, str );
end
clipboard(‘copy’,str);
end

However, when my data is as below (cellarray) in which the contents of all columns in the cellarray are numbers, but the very last column ‘info" are charss

I get this when pasted into excel (notice the last column is incorrect and actually split into 2)

this is how I created my dummy data:

uit=app.UITableSummary;
dx=9437;
z1=5795;
z2=5775.2
dz=z2-z1
tilt_mrads=-2.1084
A={dx,z1,z2,dz,tilt_mrads,’info’};
d=uit.Data; % Get Current table data
uit.Data = [d; A]; %Add new row to current data

I kept on adding this and just manually replaced the ‘info’ value in the last column (’10’, ’12’, etc..) and the numeric values in the 2nd from last column (-2.11, -1.59, etc..)

thanks for any help uitable, sprintf MATLAB Answers — New Questions

​

I have code for a battery controler i’m making, and as part of that i want to display some live results as it continues running. The whole system runs in a repreating while loop which you can see below(have removed some code to make it legible).
while On = ‘true’

%% Initialise Display
Display = figure(‘Name’,’Controler Display’,’NumberTitle’,’off’);
tiledlayout(2,1)
ax1 = nexttile;
plot(ax1,Display_timebase,Display_Power_cmd, "DisplayName","ESS Power Command (kW)")
title ‘ESS input/output’
hold on

plot(ax1,Display_timebase,Display_Grid_Demand,"DisplayName","Grid Demand")
plot(ax1,Display_timebase,Display_Available_Power, "DisplayName","Available Power (kW)")
ylim([-300 320])
legend
xline(0)
hold off
ylabel("ESS Power in/output")
xlabel("Time")
xlim([-12 48])

ax2 = nexttile;
plot(ax2,Display_timebase,Display_ESS_SoC)
title ‘SoC over time’
ylim([0 1.05])
ylabel("ESS SoC")
xlabel("Time")
xlim([-12 48])
xline(0)

linkdata(Display,’on’)
drawnow

%%%%%%%%%%%
%% Additional Code
%%%%%%%%%%%

%% Update Display values
Display_Power_cmd = [Historical_power_cmd Navigator_Outputs.Power_cmd’];
Display_Grid_Demand = [Historical_Grid_Demand Navigator_Outputs.Grid_Demand’];
Display_Available_Power = [Historical_Available_Power Mapper_Outputs.Available_Power’];
Display_ESS_SoC = [Historical_SoC ESS_SoC_array’];

drawnow

end

I’m aware that being in a while loop suppresses figure updates, hence why i’ve included the drawnows but my figures still refuse to update. What am i doing wrong here?I have code for a battery controler i’m making, and as part of that i want to display some live results as it continues running. The whole system runs in a repreating while loop which you can see below(have removed some code to make it legible).
while On = ‘true’

%% Initialise Display
Display = figure(‘Name’,’Controler Display’,’NumberTitle’,’off’);
tiledlayout(2,1)
ax1 = nexttile;
plot(ax1,Display_timebase,Display_Power_cmd, "DisplayName","ESS Power Command (kW)")
title ‘ESS input/output’
hold on

plot(ax1,Display_timebase,Display_Grid_Demand,"DisplayName","Grid Demand")
plot(ax1,Display_timebase,Display_Available_Power, "DisplayName","Available Power (kW)")
ylim([-300 320])
legend
xline(0)
hold off
ylabel("ESS Power in/output")
xlabel("Time")
xlim([-12 48])

ax2 = nexttile;
plot(ax2,Display_timebase,Display_ESS_SoC)
title ‘SoC over time’
ylim([0 1.05])
ylabel("ESS SoC")
xlabel("Time")
xlim([-12 48])
xline(0)

linkdata(Display,’on’)
drawnow

%%%%%%%%%%%
%% Additional Code
%%%%%%%%%%%

%% Update Display values
Display_Power_cmd = [Historical_power_cmd Navigator_Outputs.Power_cmd’];
Display_Grid_Demand = [Historical_Grid_Demand Navigator_Outputs.Grid_Demand’];
Display_Available_Power = [Historical_Available_Power Mapper_Outputs.Available_Power’];
Display_ESS_SoC = [Historical_SoC ESS_SoC_array’];

drawnow

end

I’m aware that being in a while loop suppresses figure updates, hence why i’ve included the drawnows but my figures still refuse to update. What am i doing wrong here? I have code for a battery controler i’m making, and as part of that i want to display some live results as it continues running. The whole system runs in a repreating while loop which you can see below(have removed some code to make it legible).
while On = ‘true’

%% Initialise Display
Display = figure(‘Name’,’Controler Display’,’NumberTitle’,’off’);
tiledlayout(2,1)
ax1 = nexttile;
plot(ax1,Display_timebase,Display_Power_cmd, "DisplayName","ESS Power Command (kW)")
title ‘ESS input/output’
hold on

plot(ax1,Display_timebase,Display_Grid_Demand,"DisplayName","Grid Demand")
plot(ax1,Display_timebase,Display_Available_Power, "DisplayName","Available Power (kW)")
ylim([-300 320])
legend
xline(0)
hold off
ylabel("ESS Power in/output")
xlabel("Time")
xlim([-12 48])

ax2 = nexttile;
plot(ax2,Display_timebase,Display_ESS_SoC)
title ‘SoC over time’
ylim([0 1.05])
ylabel("ESS SoC")
xlabel("Time")
xlim([-12 48])
xline(0)

linkdata(Display,’on’)
drawnow

%%%%%%%%%%%
%% Additional Code
%%%%%%%%%%%

%% Update Display values
Display_Power_cmd = [Historical_power_cmd Navigator_Outputs.Power_cmd’];
Display_Grid_Demand = [Historical_Grid_Demand Navigator_Outputs.Grid_Demand’];
Display_Available_Power = [Historical_Available_Power Mapper_Outputs.Available_Power’];
Display_ESS_SoC = [Historical_SoC ESS_SoC_array’];

drawnow

end

I’m aware that being in a while loop suppresses figure updates, hence why i’ve included the drawnows but my figures still refuse to update. What am i doing wrong here? figure, linkdata MATLAB Answers — New Questions

​

I implemented one custom SAC agent, which I have to, with MATLAB deep learning automatic differentiation. However, when compared to MATLAB built-in SAC agent on a certain task with exactly the same hyperparameters, the custom SAC agent failed to complete the task while the built-in agent succeeded.
Here is the training process of the built-in agent:

This is the training progress of the custom SAC agent(alongwith loss):

Here are the codes for the custom SAC agent and training:
1.Implementation of custom SAC agent
classdef MySACAgent < rl.agent.CustomAgent
properties
%networks
actor
critic1
critic2
critic_target1
critic_target2
log_alpha%entropy weight(log transformed)

%training options
options%Agent options

%optimizers
actorOptimizer
criticOptimizer_1
criticOptimizer_2
entWgtOptimizer

%experience buffers
obsBuffer
actionBuffer
rewardBuffer
nextObsBuffer
isDoneBuffer
rlExpBuffer
bufferIdx
bufferLen

%loss to record
cLoss
aLoss
eLoss

end
properties(Access = private)
Ts
counter
numObs
numAct
end

methods
%constructor
function obj = MySACAgent(numObs,numAct,obsInfo,actInfo,hid_dim,Ts,options)
% options’ field：MaxBufferLen WarmUpSteps MiniBatchSize
% LearningFrequency EntropyLossWeight DiscountFactor
% OptimizerOptions（cell） PolicyUpdateFrequency TargetEntropy
% TargetUpdateFrequency TargetSmoothFactor
% base_seed NumGradientStepsPerUpdate

%OptimizerOptions（for actor&critic）
% (required) Call the abstract class constructor.
rng(options.base_seed);%set random seed
obj = obj@rl.agent.CustomAgent();
obj.ObservationInfo = obsInfo;
obj.ActionInfo = actInfo;

% obj.SampleTime = Ts;%explicitly assigned for simulink
obj.Ts = Ts;

%create networks
if isempty(hid_dim)
hid_dim = 256;
end
obj.actor = CreateActor(obj,numObs,numAct,hid_dim,obsInfo,actInfo);
[obj.critic1,obj.critic2,obj.critic_target1,obj.critic_target2] = CreateCritic(obj,numObs,numAct,hid_dim,obsInfo,actInfo);

obj.options = options;
assert(options.WarmUpSteps>options.MiniBatchSize,…
‘options.WarmUpSteps must not be less than options.MiniBatchSize’);

%set optimizers
obj.actorOptimizer = rlOptimizer(options.OptimizerOptions{1});
obj.criticOptimizer_1 = rlOptimizer(options.OptimizerOptions{2});
obj.criticOptimizer_2 = rlOptimizer(options.OptimizerOptions{3});
obj.entWgtOptimizer = rlOptimizer(options.OptimizerOptions{4});

obj.cLoss=0;
obj.aLoss=0;
obj.eLoss=0;

% (optional) Cache the number of observations and actions.
obj.numObs = numObs;
obj.numAct = numAct;

% (optional) Initialize buffer and counter.
resetImpl(obj);
% obj.rlExpBuffer = rlReplayMemory(obsInfo,actInfo,options.MaxBufferLen);
end

function resetImpl(obj)
% (Optional) Define how the agent is reset before training/
resetBuffer(obj);
obj.counter = 0;
obj.bufferLen=0;
obj.bufferIdx = 0;%base 0
obj.log_alpha = dlarray(log(obj.options.EntropyLossWeight));
end

function resetBuffer(obj)
% Reinitialize observation buffer. Allocate as dlarray to
% support automatic differentiation with dlfeval and
% dlgradient.

%format：CBT

obj.obsBuffer = dlarray(…
zeros(obj.numObs,obj.options.MaxBufferLen),’CB’);

% Reinitialize action buffer with valid actions.
obj.actionBuffer = dlarray(…
zeros(obj.numAct,obj.options.MaxBufferLen),’CB’);

% Reinitialize reward buffer.
obj.rewardBuffer = dlarray(zeros(1,obj.options.MaxBufferLen),’CB’);

% Reinitialize nextState buffer.
obj.nextObsBuffer = dlarray(…
zeros(obj.numObs,obj.options.MaxBufferLen),’CB’);

% Reinitialize mask buffer.
obj.isDoneBuffer = dlarray(zeros(1,obj.options.MaxBufferLen),’CB’);

end

%Create networks
%Actor
function actor = CreateActor(obj,numObs,numAct,hid_dim,obsInfo,actInfo)
% Create the actor network layers.
commonPath = [
featureInputLayer(numObs,Name="obsInLyr")
fullyConnectedLayer(hid_dim)
layerNormalizationLayer
reluLayer
fullyConnectedLayer(hid_dim)
layerNormalizationLayer
reluLayer(Name="comPathOutLyr")
];
meanPath = [
fullyConnectedLayer(numAct,Name="meanOutLyr")
];
stdPath = [
fullyConnectedLayer(numAct,Name="stdInLyr")
softplusLayer(Name="stdOutLyr")
];
% Connect the layers.
actorNetwork = layerGraph(commonPath);
actorNetwork = addLayers(actorNetwork,meanPath);
actorNetwork = addLayers(actorNetwork,stdPath);
actorNetwork = connectLayers(actorNetwork,"comPathOutLyr","meanOutLyr/in");
actorNetwork = connectLayers(actorNetwork,"comPathOutLyr","stdInLyr/in");

actordlnet = dlnetwork(actorNetwork);
actor = initialize(actordlnet);
end

%Critic
function [critic1,critic2,critic_target1,critic_target2] = CreateCritic(obj,numObs,numAct,hid_dim,obsInfo,actInfo)
% Define the network layers.
criticNet = [
featureInputLayer(numObs+numAct,Name="obsInLyr")%input:[obs act]
fullyConnectedLayer(hid_dim)
layerNormalizationLayer
reluLayer
fullyConnectedLayer(hid_dim)
layerNormalizationLayer
reluLayer
fullyConnectedLayer(1,Name="QValueOutLyr")
];

% Connect the layers.
criticNet = layerGraph(criticNet);
criticDLnet = dlnetwork(criticNet,’Initialize’,false);
critic1 = initialize(criticDLnet);
critic2 = initialize(criticDLnet);%c1 and c2 different initilization
critic_target1 = initialize(criticDLnet);
critic_target1.Learnables = critic1.Learnables;
critic_target1.State = critic1.State;
critic_target2 = initialize(criticDLnet);
critic_target2.Learnables = critic2.Learnables;
critic_target2.State = critic2.State;

end

function logP = logProbBoundedAction(obj,boundedAction,mu,sigma)
%used to calculate log probability for tanh(gaussian)
%validated, nothing wrong with this function
eps=1e-10;
logP = sum(log(1/sqrt(2*pi)./sigma.*exp(-0.5*(0.5*…
log((1+boundedAction+eps)./(1-boundedAction+eps))-mu).^2./sigma.^2).*1./(1-boundedAction.^2+eps)),1);
end

%loss functions
function [vLoss_1, vLoss_2, criticGrad_1, criticGrad_2] = criticLoss(obj,batchExperiences,c1,c2)
batchObs = batchExperiences{1};
batchAction = batchExperiences{2};
batchReward = batchExperiences{3};
batchNextObs = batchExperiences{4};
batchIsDone = batchExperiences{5};
batchSize = size(batchObs,2);

gamma = obj.options.DiscountFactor;
y = dlarray(zeros(1,batchSize));%CB(C=1)
y = y + batchReward;
actionNext = getActionWithExploration_dlarray(obj,batchNextObs);%CB
actionNext = actionNext{1};
Qt1=predict(obj.critic_target1,cat(1,batchNextObs,actionNext));%CB(C=1)
Qt2=predict(obj.critic_target2,cat(1,batchNextObs,actionNext));%CB(C=1)

[mu,sigma] = predict(obj.actor,batchNextObs);%CB:numAct*batch
next_action = tanh(mu + sigma.*randn(size(sigma)));
logP = logProbBoundedAction(obj,next_action,mu,sigma);
y = y + (1 – batchIsDone).*(gamma*(min(cat(1,Qt1,Qt2),[],1) – exp(obj.log_alpha)*logP));

critic_input = cat(1,batchObs,batchAction);
Q1 = forward(c1,critic_input);
Q2 = forward(c2,critic_input);
vLoss_1 = 1/2*mean((y – Q1).^2,’all’);
vLoss_2 = 1/2*mean((y – Q2).^2,’all’);
criticGrad_1 = dlgradient(vLoss_1,c1.Learnables);
criticGrad_2 = dlgradient(vLoss_2,c2.Learnables);
end

function [aLoss,actorGrad] = actorLoss(obj,batchExperiences,actor)
batchObs = batchExperiences{1};
batchSize = size(batchObs,2);
[mu,sigma] = forward(actor,batchObs);%CB:numAct*batch
curr_action = tanh(mu + sigma.*randn(size(sigma)));%reparameterization
critic_input = cat(1,batchObs,curr_action);
Q1=forward(obj.critic1,critic_input);%CB(C=1)
Q2=forward(obj.critic2,critic_input);%CB(C=1)

logP = logProbBoundedAction(obj,curr_action,mu,sigma);
aLoss = mean(-min(cat(1,Q1,Q2),[],1) + exp(obj.log_alpha) * logP,’all’);

actorGrad= dlgradient(aLoss,actor.Learnables);
end

function [eLoss,entGrad] = entropyLoss(obj,batchExperiences,logAlpha)
batchObs = batchExperiences{1};
[mu,sigma] = predict(obj.actor,batchObs);%CB:numAct*batch
curr_action = tanh(mu + sigma.*randn(size(sigma)));
ent = mean(-logProbBoundedAction(obj,curr_action,mu,sigma));
eLoss = exp(logAlpha) * (ent – obj.options.TargetEntropy);
entGrad = dlgradient(eLoss,logAlpha);
end

end
methods(Access=protected)
%return SampleTime
function ts = getSampleTime_(obj)
ts = obj.Ts;
end
%get action without exploration
function action = getActionImpl(obj,obs)
%obs:dlarray CB
if ~isa(obs,’dlarray’)
if isa(obs,’cell’)
obs = dlarray(obs{1},’CB’);
else
obs = dlarray(obs,’CB’);
end
end
[mu,~] = predict(obj.actor,obs);
mu = extractdata(mu);
action = {tanh(mu)};
end

%get action with exploration
function action = getActionWithExplorationImpl(obj,obs)
%obs:dlarray CT
if ~isa(obs,’dlarray’) || size(obs,1)~=obj.numObs
obs = dlarray(randn(obj.numObs,1),’CB’);
end
[mu,sigma] = predict(obj.actor,obs);
mu = extractdata(mu);
sigma = extractdata(sigma);
action = {tanh(mu + sigma .* randn(size(sigma)))};
end

function action = getActionWithExploration_dlarray(obj,obs)
[mu,sigma] = predict(obj.actor,obs);
action = {tanh(mu + sigma .* randn(size(sigma)))};
end

%learning
function action = learnImpl(obj,Experience)
% Extract data from experience.

obs = Experience{1};
action = Experience{2};
reward = Experience{3};
nextObs = Experience{4};
isDone = logical(Experience{5});

obj.obsBuffer(:,obj.bufferIdx+1,:) = obs{1};
obj.actionBuffer(:,obj.bufferIdx+1,:) = action{1};
obj.rewardBuffer(:,obj.bufferIdx+1) = reward;
obj.nextObsBuffer(:,obj.bufferIdx+1,:) = nextObs{1};
obj.isDoneBuffer(:,obj.bufferIdx+1) = isDone;
obj.bufferLen = max(obj.bufferLen,obj.bufferIdx+1);
obj.bufferIdx = mod(obj.bufferIdx+1,obj.options.MaxBufferLen);

if obj.bufferLen>=max(obj.options.WarmUpSteps,obj.options.MiniBatchSize)
obj.counter = obj.counter + 1;
if (obj.options.LearningFrequency==-1 && isDone) || …
(obj.options.LearningFrequency>0 && mod(obj.counter,obj.options.LearningFrequency)==0)
for gstep = 1:obj.options.NumGradientStepsPerUpdate
%sample batch
batchSize = obj.options.MiniBatchSize;
batchInd = randperm(obj.bufferLen,batchSize);

batchExperience = {
obj.obsBuffer(:,batchInd,:),…
obj.actionBuffer(:,batchInd,:),…
obj.rewardBuffer(:,batchInd),…
obj.nextObsBuffer(:,batchInd,:),…
obj.isDoneBuffer(:,batchInd)
};

%update the parameters of each critic
[cLoss1,cLoss2,criticGrad_1,criticGrad_2] = dlfeval(@(x,c1,c2)obj.criticLoss(x,c1,c2),batchExperience,obj.critic1,obj.critic2);
obj.cLoss = min(extractdata(cLoss1),extractdata(cLoss2));
[obj.critic1.Learnables.Value,obj.criticOptimizer_1] = update(obj.criticOptimizer_1,obj.critic1.Learnables.Value,criticGrad_1.Value);
[obj.critic2.Learnables.Value,obj.criticOptimizer_2] = update(obj.criticOptimizer_2,obj.critic2.Learnables.Value,criticGrad_2.Value);

if (mod(obj.counter,obj.options.PolicyUpdateFrequency)==0 && obj.options.LearningFrequency==-1) ||…
(mod(obj.counter,obj.options.LearningFrequency * obj.options.PolicyUpdateFrequency)==0 …
&& obj.options.LearningFrequency>0)
%update the parameters of actor
[aloss,actorGrad] = dlfeval(…
@(x,actor)obj.actorLoss(x,actor),…
batchExperience,obj.actor);
obj.aLoss = extractdata(aloss);
[obj.actor.Learnables.Value,obj.actorOptimizer] = update(obj.actorOptimizer,obj.actor.Learnables.Value,actorGrad.Value);

%update the entropy weight
[eloss,entGrad] = dlfeval(@(x,alpha)obj.entropyLoss(x,alpha),batchExperience,obj.log_alpha);
obj.eLoss = extractdata(eloss);
% disp(obj.alpha)
[obj.log_alpha,obj.entWgtOptimizer] = update(obj.entWgtOptimizer,{obj.log_alpha},{entGrad});
obj.log_alpha = obj.log_alpha{1};
end

%update critic targets
%1
critic1_params = obj.critic1.Learnables.Value;%cell array network params
critic_target1_params = obj.critic_target1.Learnables.Value;
for i=1:size(critic1_params,1)
obj.critic_target1.Learnables.Value{i} = obj.options.TargetSmoothFactor * critic1_params{i}…
+ (1 – obj.options.TargetSmoothFactor) * critic_target1_params{i};
end
%2
critic2_params = obj.critic2.Learnables.Value;%cell array network params
critic_target2_params = obj.critic_target2.Learnables.Value;
for i=1:size(critic2_params,1)
obj.critic_target2.Learnables.Value{i} = obj.options.TargetSmoothFactor * critic2_params{i}…
+ (1 – obj.options.TargetSmoothFactor) * critic_target2_params{i};
end
% end
end
end

end
action = getActionWithExplorationImpl(obj,nextObs{1});
end
end
end

2.Configuration of ‘options’ property(same as those used for the built-in SAC agent)
options.MaxBufferLen = 1e4;
options.WarmUpSteps = 1000;
options.MiniBatchSize = 256;
options.LearningFrequency = -1;%when -1: train after each episode
options.EntropyLossWeight = 1;
options.DiscountFactor = 0.99;
options.PolicyUpdateFrequency = 1;
options.TargetEntropy = -2;
options.TargetUpdateFrequency = 1;
options.TargetSmoothFactor = 1e-3;
options.NumGradientStepsPerUpdate = 10;
%optimizerOptions: actor critic1 critic2 entWgt(alpha)
%encoder decoder
options.OptimizerOptions = {
rlOptimizerOptions("Algorithm","adam","GradientThreshold",1,’LearnRate’,1e-3),…
rlOptimizerOptions("Algorithm","adam","GradientThreshold",1,’LearnRate’,1e-3),…
rlOptimizerOptions("Algorithm","adam","GradientThreshold",1,’LearnRate’,1e-3),…
rlOptimizerOptions("Algorithm","adam",’LearnRate’,3e-4),…
rlOptimizerOptions("Algorithm","adam","GradientThreshold",1,’LearnRate’,1e-3),…
rlOptimizerOptions("Algorithm","adam","GradientThreshold",1,’LearnRate’,1e-3)};
options.base_seed=940;
3.training

clc;
clear;
close all;
run(‘init_car_params.m’);
%create RL env
numObs = 4; % vx vy r beta_user
numAct = 2; % st_angle_ref rw_omega_ref

obsInfo = rlNumericSpec([numObs 1]);

actInfo = rlNumericSpec([numAct 1]);
actInfo.LowerLimit = -1;
actInfo.UpperLimit = 1;

mdl = "prius_sm_model";
blk = mdl + "/RL Agent";
env = rlSimulinkEnv(mdl,blk,obsInfo,actInfo);

params=struct(‘rw_radius’,rw_radius,’a’,a,’b’,b,’init_vx’,init_vx,’init_yaw_rate’,init_yaw_rate);
env.ResetFcn = @(in) PriusResetFcn(in,params,mdl);

Ts = 1/10;
Tf = 5;

%create actor
rnd_seed=940;
algorithm = ‘MySAC’;
switch algorithm
case ‘SAC’
agent = createNetworks(rnd_seed,numObs,numAct,obsInfo,actInfo,Ts);
case ‘MySAC’
hid_dim = 256;
options=getDWMLAgentOptions();
agent = MySACAgent(numObs,numAct,obsInfo,actInfo,hid_dim,Ts,options);
end

%%
%train agent
close all
maxEpisodes = 6000;
maxSteps = floor(Tf/Ts);
useParallel = false;
run_idx=9;
saveAgentDir = [‘savedAgents/’,algorithm,’/’,num2str(run_idx)];
switch algorithm
case ‘SAC’
trainOpts = rlTrainingOptions(…
MaxEpisodes=maxEpisodes, …
MaxStepsPerEpisode=maxSteps, …
ScoreAveragingWindowLength=100, …
Plots="training-progress", …
StopTrainingCriteria="AverageReward", …
UseParallel=useParallel,…
SaveAgentCriteria=’EpisodeReward’,…
SaveAgentValue=35,…
SaveAgentDirectory=saveAgentDir);
% SaveAgentCriteria=’EpisodeFrequency’,…
% SaveAgentValue=1,…
case ‘MySAC’
trainOpts = rlTrainingOptions(…
MaxEpisodes=maxEpisodes, …
MaxStepsPerEpisode=maxSteps, …
ScoreAveragingWindowLength=100, …
Plots="training-progress", …
StopTrainingCriteria="AverageReward", …
UseParallel=useParallel,…
SaveAgentCriteria=’EpisodeReward’,…
SaveAgentValue=35,…
SaveAgentDirectory=saveAgentDir);
end

set_param(mdl,"FastRestart","off");%for random initialization
if trainOpts.UseParallel
% Disable visualization in Simscape Mechanics Explorer
set_param(mdl, SimMechanicsOpenEditorOnUpdate="off");
save_system(mdl);
else
% Enable visualization in Simscape Mechanics Explorer
set_param(mdl, SimMechanicsOpenEditorOnUpdate="on");
end

%load training data
monitor = trainingProgressMonitor();
logger = rlDataLogger(monitor);
logger.EpisodeFinishedFcn = @myEpisodeLoggingFcn;

doTraining = true;
if doTraining
trainResult = train(agent,env,trainOpts,Logger=logger);
end

% %logger callback used for MySACAgent
function dataToLog = myEpisodeLoggingFcn(data)
dataToLog.criticLoss = data.Agent.cLoss;
dataToLog.actorLoss = data.Agent.aLoss;
dataToLog.entLoss = data.Agent.eLoss;
% dataToLog.denoiseLoss = data.Agent.dnLoss;
end
In the simulink environment used, action output by the Agent block(in [-1,1]) is denormalized and fed into the environment.

I think possible causes of the problem include:
1.Wrong implementation of critic loss. As shown in the training progress, critic loss seemed to diverge. It’s hardly caused by hyperparameters(batch size or learning rate or target update frequency) because they worked well for the built-in agent. So it is more likely the critic loss is wrong.
2.Wrong implementation of replay buffer. I implemented the replay buffer as a circular queue, where I sampled uniformly to get batch training data. From the comparison of the training progress shown above, the custom SAC agent did explore states with high reward(around 30) but failed to exploit them, So I guess there is still problem with my replay buffer.
3.Gradient flow was broken.The learning is done with the help of MATLAB deep learning automatic differentiation. Perhaps some of my implementation violates the computational rule of automatic differentiation, which broke the gradient flow during forward computation or backpropagation and led to wrong result.
4.Gradient step(update frequency). In current implementation, NumGradientStepsPerUpdate gradient steps are executed after each episode. During each gradient step, cirtic(s) and actor, alongwith entropy weight, is updated once. I am not sure whether the current implementation of gradient step has got the update frequency right.
5.Also could be normalization problem, but I am not so sure.
I plan to debug 3 first.
Please read the code and help find potential causes of the gap between the custom SAC agent and the built-in one.
Finally, I am actually trying to extend SAC algorithm to a more complex framework. I didn’t choose to inherit the built-in SAC agent(rlSACAgent), would it be recommended to do my development by doing so?I implemented one custom SAC agent, which I have to, with MATLAB deep learning automatic differentiation. However, when compared to MATLAB built-in SAC agent on a certain task with exactly the same hyperparameters, the custom SAC agent failed to complete the task while the built-in agent succeeded.
Here is the training process of the built-in agent:

This is the training progress of the custom SAC agent(alongwith loss):

Here are the codes for the custom SAC agent and training:
1.Implementation of custom SAC agent
classdef MySACAgent < rl.agent.CustomAgent
properties
%networks
actor
critic1
critic2
critic_target1
critic_target2
log_alpha%entropy weight(log transformed)

%training options
options%Agent options

%optimizers
actorOptimizer
criticOptimizer_1
criticOptimizer_2
entWgtOptimizer

%experience buffers
obsBuffer
actionBuffer
rewardBuffer
nextObsBuffer
isDoneBuffer
rlExpBuffer
bufferIdx
bufferLen

%loss to record
cLoss
aLoss
eLoss

end
properties(Access = private)
Ts
counter
numObs
numAct
end

methods
%constructor
function obj = MySACAgent(numObs,numAct,obsInfo,actInfo,hid_dim,Ts,options)
% options’ field：MaxBufferLen WarmUpSteps MiniBatchSize
% LearningFrequency EntropyLossWeight DiscountFactor
% OptimizerOptions（cell） PolicyUpdateFrequency TargetEntropy
% TargetUpdateFrequency TargetSmoothFactor
% base_seed NumGradientStepsPerUpdate

%OptimizerOptions（for actor&critic）
% (required) Call the abstract class constructor.
rng(options.base_seed);%set random seed
obj = obj@rl.agent.CustomAgent();
obj.ObservationInfo = obsInfo;
obj.ActionInfo = actInfo;

% obj.SampleTime = Ts;%explicitly assigned for simulink
obj.Ts = Ts;

%create networks
if isempty(hid_dim)
hid_dim = 256;
end
obj.actor = CreateActor(obj,numObs,numAct,hid_dim,obsInfo,actInfo);
[obj.critic1,obj.critic2,obj.critic_target1,obj.critic_target2] = CreateCritic(obj,numObs,numAct,hid_dim,obsInfo,actInfo);

obj.options = options;
assert(options.WarmUpSteps>options.MiniBatchSize,…
‘options.WarmUpSteps must not be less than options.MiniBatchSize’);

%set optimizers
obj.actorOptimizer = rlOptimizer(options.OptimizerOptions{1});
obj.criticOptimizer_1 = rlOptimizer(options.OptimizerOptions{2});
obj.criticOptimizer_2 = rlOptimizer(options.OptimizerOptions{3});
obj.entWgtOptimizer = rlOptimizer(options.OptimizerOptions{4});

obj.cLoss=0;
obj.aLoss=0;
obj.eLoss=0;

% (optional) Cache the number of observations and actions.
obj.numObs = numObs;
obj.numAct = numAct;

% (optional) Initialize buffer and counter.
resetImpl(obj);
% obj.rlExpBuffer = rlReplayMemory(obsInfo,actInfo,options.MaxBufferLen);
end

function resetImpl(obj)
% (Optional) Define how the agent is reset before training/
resetBuffer(obj);
obj.counter = 0;
obj.bufferLen=0;
obj.bufferIdx = 0;%base 0
obj.log_alpha = dlarray(log(obj.options.EntropyLossWeight));
end

function resetBuffer(obj)
% Reinitialize observation buffer. Allocate as dlarray to
% support automatic differentiation with dlfeval and
% dlgradient.

%format：CBT

obj.obsBuffer = dlarray(…
zeros(obj.numObs,obj.options.MaxBufferLen),’CB’);

% Reinitialize action buffer with valid actions.
obj.actionBuffer = dlarray(…
zeros(obj.numAct,obj.options.MaxBufferLen),’CB’);

% Reinitialize reward buffer.
obj.rewardBuffer = dlarray(zeros(1,obj.options.MaxBufferLen),’CB’);

% Reinitialize nextState buffer.
obj.nextObsBuffer = dlarray(…
zeros(obj.numObs,obj.options.MaxBufferLen),’CB’);

% Reinitialize mask buffer.
obj.isDoneBuffer = dlarray(zeros(1,obj.options.MaxBufferLen),’CB’);

end

%Create networks
%Actor
function actor = CreateActor(obj,numObs,numAct,hid_dim,obsInfo,actInfo)
% Create the actor network layers.
commonPath = [
featureInputLayer(numObs,Name="obsInLyr")
fullyConnectedLayer(hid_dim)
layerNormalizationLayer
reluLayer
fullyConnectedLayer(hid_dim)
layerNormalizationLayer
reluLayer(Name="comPathOutLyr")
];
meanPath = [
fullyConnectedLayer(numAct,Name="meanOutLyr")
];
stdPath = [
fullyConnectedLayer(numAct,Name="stdInLyr")
softplusLayer(Name="stdOutLyr")
];
% Connect the layers.
actorNetwork = layerGraph(commonPath);
actorNetwork = addLayers(actorNetwork,meanPath);
actorNetwork = addLayers(actorNetwork,stdPath);
actorNetwork = connectLayers(actorNetwork,"comPathOutLyr","meanOutLyr/in");
actorNetwork = connectLayers(actorNetwork,"comPathOutLyr","stdInLyr/in");

actordlnet = dlnetwork(actorNetwork);
actor = initialize(actordlnet);
end

%Critic
function [critic1,critic2,critic_target1,critic_target2] = CreateCritic(obj,numObs,numAct,hid_dim,obsInfo,actInfo)
% Define the network layers.
criticNet = [
featureInputLayer(numObs+numAct,Name="obsInLyr")%input:[obs act]
fullyConnectedLayer(hid_dim)
layerNormalizationLayer
reluLayer
fullyConnectedLayer(hid_dim)
layerNormalizationLayer
reluLayer
fullyConnectedLayer(1,Name="QValueOutLyr")
];

% Connect the layers.
criticNet = layerGraph(criticNet);
criticDLnet = dlnetwork(criticNet,’Initialize’,false);
critic1 = initialize(criticDLnet);
critic2 = initialize(criticDLnet);%c1 and c2 different initilization
critic_target1 = initialize(criticDLnet);
critic_target1.Learnables = critic1.Learnables;
critic_target1.State = critic1.State;
critic_target2 = initialize(criticDLnet);
critic_target2.Learnables = critic2.Learnables;
critic_target2.State = critic2.State;

end

function logP = logProbBoundedAction(obj,boundedAction,mu,sigma)
%used to calculate log probability for tanh(gaussian)
%validated, nothing wrong with this function
eps=1e-10;
logP = sum(log(1/sqrt(2*pi)./sigma.*exp(-0.5*(0.5*…
log((1+boundedAction+eps)./(1-boundedAction+eps))-mu).^2./sigma.^2).*1./(1-boundedAction.^2+eps)),1);
end

%loss functions
function [vLoss_1, vLoss_2, criticGrad_1, criticGrad_2] = criticLoss(obj,batchExperiences,c1,c2)
batchObs = batchExperiences{1};
batchAction = batchExperiences{2};
batchReward = batchExperiences{3};
batchNextObs = batchExperiences{4};
batchIsDone = batchExperiences{5};
batchSize = size(batchObs,2);

gamma = obj.options.DiscountFactor;
y = dlarray(zeros(1,batchSize));%CB(C=1)
y = y + batchReward;
actionNext = getActionWithExploration_dlarray(obj,batchNextObs);%CB
actionNext = actionNext{1};
Qt1=predict(obj.critic_target1,cat(1,batchNextObs,actionNext));%CB(C=1)
Qt2=predict(obj.critic_target2,cat(1,batchNextObs,actionNext));%CB(C=1)

[mu,sigma] = predict(obj.actor,batchNextObs);%CB:numAct*batch
next_action = tanh(mu + sigma.*randn(size(sigma)));
logP = logProbBoundedAction(obj,next_action,mu,sigma);
y = y + (1 – batchIsDone).*(gamma*(min(cat(1,Qt1,Qt2),[],1) – exp(obj.log_alpha)*logP));

critic_input = cat(1,batchObs,batchAction);
Q1 = forward(c1,critic_input);
Q2 = forward(c2,critic_input);
vLoss_1 = 1/2*mean((y – Q1).^2,’all’);
vLoss_2 = 1/2*mean((y – Q2).^2,’all’);
criticGrad_1 = dlgradient(vLoss_1,c1.Learnables);
criticGrad_2 = dlgradient(vLoss_2,c2.Learnables);
end

function [aLoss,actorGrad] = actorLoss(obj,batchExperiences,actor)
batchObs = batchExperiences{1};
batchSize = size(batchObs,2);
[mu,sigma] = forward(actor,batchObs);%CB:numAct*batch
curr_action = tanh(mu + sigma.*randn(size(sigma)));%reparameterization
critic_input = cat(1,batchObs,curr_action);
Q1=forward(obj.critic1,critic_input);%CB(C=1)
Q2=forward(obj.critic2,critic_input);%CB(C=1)

logP = logProbBoundedAction(obj,curr_action,mu,sigma);
aLoss = mean(-min(cat(1,Q1,Q2),[],1) + exp(obj.log_alpha) * logP,’all’);

actorGrad= dlgradient(aLoss,actor.Learnables);
end

function [eLoss,entGrad] = entropyLoss(obj,batchExperiences,logAlpha)
batchObs = batchExperiences{1};
[mu,sigma] = predict(obj.actor,batchObs);%CB:numAct*batch
curr_action = tanh(mu + sigma.*randn(size(sigma)));
ent = mean(-logProbBoundedAction(obj,curr_action,mu,sigma));
eLoss = exp(logAlpha) * (ent – obj.options.TargetEntropy);
entGrad = dlgradient(eLoss,logAlpha);
end

end
methods(Access=protected)
%return SampleTime
function ts = getSampleTime_(obj)
ts = obj.Ts;
end
%get action without exploration
function action = getActionImpl(obj,obs)
%obs:dlarray CB
if ~isa(obs,’dlarray’)
if isa(obs,’cell’)
obs = dlarray(obs{1},’CB’);
else
obs = dlarray(obs,’CB’);
end
end
[mu,~] = predict(obj.actor,obs);
mu = extractdata(mu);
action = {tanh(mu)};
end

%get action with exploration
function action = getActionWithExplorationImpl(obj,obs)
%obs:dlarray CT
if ~isa(obs,’dlarray’) || size(obs,1)~=obj.numObs
obs = dlarray(randn(obj.numObs,1),’CB’);
end
[mu,sigma] = predict(obj.actor,obs);
mu = extractdata(mu);
sigma = extractdata(sigma);
action = {tanh(mu + sigma .* randn(size(sigma)))};
end

function action = getActionWithExploration_dlarray(obj,obs)
[mu,sigma] = predict(obj.actor,obs);
action = {tanh(mu + sigma .* randn(size(sigma)))};
end

%learning
function action = learnImpl(obj,Experience)
% Extract data from experience.

obs = Experience{1};
action = Experience{2};
reward = Experience{3};
nextObs = Experience{4};
isDone = logical(Experience{5});

obj.obsBuffer(:,obj.bufferIdx+1,:) = obs{1};
obj.actionBuffer(:,obj.bufferIdx+1,:) = action{1};
obj.rewardBuffer(:,obj.bufferIdx+1) = reward;
obj.nextObsBuffer(:,obj.bufferIdx+1,:) = nextObs{1};
obj.isDoneBuffer(:,obj.bufferIdx+1) = isDone;
obj.bufferLen = max(obj.bufferLen,obj.bufferIdx+1);
obj.bufferIdx = mod(obj.bufferIdx+1,obj.options.MaxBufferLen);

if obj.bufferLen>=max(obj.options.WarmUpSteps,obj.options.MiniBatchSize)
obj.counter = obj.counter + 1;
if (obj.options.LearningFrequency==-1 && isDone) || …
(obj.options.LearningFrequency>0 && mod(obj.counter,obj.options.LearningFrequency)==0)
for gstep = 1:obj.options.NumGradientStepsPerUpdate
%sample batch
batchSize = obj.options.MiniBatchSize;
batchInd = randperm(obj.bufferLen,batchSize);

batchExperience = {
obj.obsBuffer(:,batchInd,:),…
obj.actionBuffer(:,batchInd,:),…
obj.rewardBuffer(:,batchInd),…
obj.nextObsBuffer(:,batchInd,:),…
obj.isDoneBuffer(:,batchInd)
};

%update the parameters of each critic
[cLoss1,cLoss2,criticGrad_1,criticGrad_2] = dlfeval(@(x,c1,c2)obj.criticLoss(x,c1,c2),batchExperience,obj.critic1,obj.critic2);
obj.cLoss = min(extractdata(cLoss1),extractdata(cLoss2));
[obj.critic1.Learnables.Value,obj.criticOptimizer_1] = update(obj.criticOptimizer_1,obj.critic1.Learnables.Value,criticGrad_1.Value);
[obj.critic2.Learnables.Value,obj.criticOptimizer_2] = update(obj.criticOptimizer_2,obj.critic2.Learnables.Value,criticGrad_2.Value);

if (mod(obj.counter,obj.options.PolicyUpdateFrequency)==0 && obj.options.LearningFrequency==-1) ||…
(mod(obj.counter,obj.options.LearningFrequency * obj.options.PolicyUpdateFrequency)==0 …
&& obj.options.LearningFrequency>0)
%update the parameters of actor
[aloss,actorGrad] = dlfeval(…
@(x,actor)obj.actorLoss(x,actor),…
batchExperience,obj.actor);
obj.aLoss = extractdata(aloss);
[obj.actor.Learnables.Value,obj.actorOptimizer] = update(obj.actorOptimizer,obj.actor.Learnables.Value,actorGrad.Value);

%update the entropy weight
[eloss,entGrad] = dlfeval(@(x,alpha)obj.entropyLoss(x,alpha),batchExperience,obj.log_alpha);
obj.eLoss = extractdata(eloss);
% disp(obj.alpha)
[obj.log_alpha,obj.entWgtOptimizer] = update(obj.entWgtOptimizer,{obj.log_alpha},{entGrad});
obj.log_alpha = obj.log_alpha{1};
end

%update critic targets
%1
critic1_params = obj.critic1.Learnables.Value;%cell array network params
critic_target1_params = obj.critic_target1.Learnables.Value;
for i=1:size(critic1_params,1)
obj.critic_target1.Learnables.Value{i} = obj.options.TargetSmoothFactor * critic1_params{i}…
+ (1 – obj.options.TargetSmoothFactor) * critic_target1_params{i};
end
%2
critic2_params = obj.critic2.Learnables.Value;%cell array network params
critic_target2_params = obj.critic_target2.Learnables.Value;
for i=1:size(critic2_params,1)
obj.critic_target2.Learnables.Value{i} = obj.options.TargetSmoothFactor * critic2_params{i}…
+ (1 – obj.options.TargetSmoothFactor) * critic_target2_params{i};
end
% end
end
end

end
action = getActionWithExplorationImpl(obj,nextObs{1});
end
end
end

2.Configuration of ‘options’ property(same as those used for the built-in SAC agent)
options.MaxBufferLen = 1e4;
options.WarmUpSteps = 1000;
options.MiniBatchSize = 256;
options.LearningFrequency = -1;%when -1: train after each episode
options.EntropyLossWeight = 1;
options.DiscountFactor = 0.99;
options.PolicyUpdateFrequency = 1;
options.TargetEntropy = -2;
options.TargetUpdateFrequency = 1;
options.TargetSmoothFactor = 1e-3;
options.NumGradientStepsPerUpdate = 10;
%optimizerOptions: actor critic1 critic2 entWgt(alpha)
%encoder decoder
options.OptimizerOptions = {
rlOptimizerOptions("Algorithm","adam","GradientThreshold",1,’LearnRate’,1e-3),…
rlOptimizerOptions("Algorithm","adam","GradientThreshold",1,’LearnRate’,1e-3),…
rlOptimizerOptions("Algorithm","adam","GradientThreshold",1,’LearnRate’,1e-3),…
rlOptimizerOptions("Algorithm","adam",’LearnRate’,3e-4),…
rlOptimizerOptions("Algorithm","adam","GradientThreshold",1,’LearnRate’,1e-3),…
rlOptimizerOptions("Algorithm","adam","GradientThreshold",1,’LearnRate’,1e-3)};
options.base_seed=940;
3.training

clc;
clear;
close all;
run(‘init_car_params.m’);
%create RL env
numObs = 4; % vx vy r beta_user
numAct = 2; % st_angle_ref rw_omega_ref

obsInfo = rlNumericSpec([numObs 1]);

actInfo = rlNumericSpec([numAct 1]);
actInfo.LowerLimit = -1;
actInfo.UpperLimit = 1;

mdl = "prius_sm_model";
blk = mdl + "/RL Agent";
env = rlSimulinkEnv(mdl,blk,obsInfo,actInfo);

params=struct(‘rw_radius’,rw_radius,’a’,a,’b’,b,’init_vx’,init_vx,’init_yaw_rate’,init_yaw_rate);
env.ResetFcn = @(in) PriusResetFcn(in,params,mdl);

Ts = 1/10;
Tf = 5;

%create actor
rnd_seed=940;
algorithm = ‘MySAC’;
switch algorithm
case ‘SAC’
agent = createNetworks(rnd_seed,numObs,numAct,obsInfo,actInfo,Ts);
case ‘MySAC’
hid_dim = 256;
options=getDWMLAgentOptions();
agent = MySACAgent(numObs,numAct,obsInfo,actInfo,hid_dim,Ts,options);
end

%%
%train agent
close all
maxEpisodes = 6000;
maxSteps = floor(Tf/Ts);
useParallel = false;
run_idx=9;
saveAgentDir = [‘savedAgents/’,algorithm,’/’,num2str(run_idx)];
switch algorithm
case ‘SAC’
trainOpts = rlTrainingOptions(…
MaxEpisodes=maxEpisodes, …
MaxStepsPerEpisode=maxSteps, …
ScoreAveragingWindowLength=100, …
Plots="training-progress", …
StopTrainingCriteria="AverageReward", …
UseParallel=useParallel,…
SaveAgentCriteria=’EpisodeReward’,…
SaveAgentValue=35,…
SaveAgentDirectory=saveAgentDir);
% SaveAgentCriteria=’EpisodeFrequency’,…
% SaveAgentValue=1,…
case ‘MySAC’
trainOpts = rlTrainingOptions(…
MaxEpisodes=maxEpisodes, …
MaxStepsPerEpisode=maxSteps, …
ScoreAveragingWindowLength=100, …
Plots="training-progress", …
StopTrainingCriteria="AverageReward", …
UseParallel=useParallel,…
SaveAgentCriteria=’EpisodeReward’,…
SaveAgentValue=35,…
SaveAgentDirectory=saveAgentDir);
end

set_param(mdl,"FastRestart","off");%for random initialization
if trainOpts.UseParallel
% Disable visualization in Simscape Mechanics Explorer
set_param(mdl, SimMechanicsOpenEditorOnUpdate="off");
save_system(mdl);
else
% Enable visualization in Simscape Mechanics Explorer
set_param(mdl, SimMechanicsOpenEditorOnUpdate="on");
end

%load training data
monitor = trainingProgressMonitor();
logger = rlDataLogger(monitor);
logger.EpisodeFinishedFcn = @myEpisodeLoggingFcn;

doTraining = true;
if doTraining
trainResult = train(agent,env,trainOpts,Logger=logger);
end

% %logger callback used for MySACAgent
function dataToLog = myEpisodeLoggingFcn(data)
dataToLog.criticLoss = data.Agent.cLoss;
dataToLog.actorLoss = data.Agent.aLoss;
dataToLog.entLoss = data.Agent.eLoss;
% dataToLog.denoiseLoss = data.Agent.dnLoss;
end
In the simulink environment used, action output by the Agent block(in [-1,1]) is denormalized and fed into the environment.

I think possible causes of the problem include:
1.Wrong implementation of critic loss. As shown in the training progress, critic loss seemed to diverge. It’s hardly caused by hyperparameters(batch size or learning rate or target update frequency) because they worked well for the built-in agent. So it is more likely the critic loss is wrong.
2.Wrong implementation of replay buffer. I implemented the replay buffer as a circular queue, where I sampled uniformly to get batch training data. From the comparison of the training progress shown above, the custom SAC agent did explore states with high reward(around 30) but failed to exploit them, So I guess there is still problem with my replay buffer.
3.Gradient flow was broken.The learning is done with the help of MATLAB deep learning automatic differentiation. Perhaps some of my implementation violates the computational rule of automatic differentiation, which broke the gradient flow during forward computation or backpropagation and led to wrong result.
4.Gradient step(update frequency). In current implementation, NumGradientStepsPerUpdate gradient steps are executed after each episode. During each gradient step, cirtic(s) and actor, alongwith entropy weight, is updated once. I am not sure whether the current implementation of gradient step has got the update frequency right.
5.Also could be normalization problem, but I am not so sure.
I plan to debug 3 first.
Please read the code and help find potential causes of the gap between the custom SAC agent and the built-in one.
Finally, I am actually trying to extend SAC algorithm to a more complex framework. I didn’t choose to inherit the built-in SAC agent(rlSACAgent), would it be recommended to do my development by doing so? I implemented one custom SAC agent, which I have to, with MATLAB deep learning automatic differentiation. However, when compared to MATLAB built-in SAC agent on a certain task with exactly the same hyperparameters, the custom SAC agent failed to complete the task while the built-in agent succeeded.
Here is the training process of the built-in agent:

This is the training progress of the custom SAC agent(alongwith loss):

Here are the codes for the custom SAC agent and training:
1.Implementation of custom SAC agent
classdef MySACAgent < rl.agent.CustomAgent
properties
%networks
actor
critic1
critic2
critic_target1
critic_target2
log_alpha%entropy weight(log transformed)

%training options
options%Agent options

%optimizers
actorOptimizer
criticOptimizer_1
criticOptimizer_2
entWgtOptimizer

%experience buffers
obsBuffer
actionBuffer
rewardBuffer
nextObsBuffer
isDoneBuffer
rlExpBuffer
bufferIdx
bufferLen

%loss to record
cLoss
aLoss
eLoss

end
properties(Access = private)
Ts
counter
numObs
numAct
end

methods
%constructor
function obj = MySACAgent(numObs,numAct,obsInfo,actInfo,hid_dim,Ts,options)
% options’ field：MaxBufferLen WarmUpSteps MiniBatchSize
% LearningFrequency EntropyLossWeight DiscountFactor
% OptimizerOptions（cell） PolicyUpdateFrequency TargetEntropy
% TargetUpdateFrequency TargetSmoothFactor
% base_seed NumGradientStepsPerUpdate

%OptimizerOptions（for actor&critic）
% (required) Call the abstract class constructor.
rng(options.base_seed);%set random seed
obj = obj@rl.agent.CustomAgent();
obj.ObservationInfo = obsInfo;
obj.ActionInfo = actInfo;

% obj.SampleTime = Ts;%explicitly assigned for simulink
obj.Ts = Ts;

%create networks
if isempty(hid_dim)
hid_dim = 256;
end
obj.actor = CreateActor(obj,numObs,numAct,hid_dim,obsInfo,actInfo);
[obj.critic1,obj.critic2,obj.critic_target1,obj.critic_target2] = CreateCritic(obj,numObs,numAct,hid_dim,obsInfo,actInfo);

obj.options = options;
assert(options.WarmUpSteps>options.MiniBatchSize,…
‘options.WarmUpSteps must not be less than options.MiniBatchSize’);

%set optimizers
obj.actorOptimizer = rlOptimizer(options.OptimizerOptions{1});
obj.criticOptimizer_1 = rlOptimizer(options.OptimizerOptions{2});
obj.criticOptimizer_2 = rlOptimizer(options.OptimizerOptions{3});
obj.entWgtOptimizer = rlOptimizer(options.OptimizerOptions{4});

obj.cLoss=0;
obj.aLoss=0;
obj.eLoss=0;

% (optional) Cache the number of observations and actions.
obj.numObs = numObs;
obj.numAct = numAct;

% (optional) Initialize buffer and counter.
resetImpl(obj);
% obj.rlExpBuffer = rlReplayMemory(obsInfo,actInfo,options.MaxBufferLen);
end

function resetImpl(obj)
% (Optional) Define how the agent is reset before training/
resetBuffer(obj);
obj.counter = 0;
obj.bufferLen=0;
obj.bufferIdx = 0;%base 0
obj.log_alpha = dlarray(log(obj.options.EntropyLossWeight));
end

function resetBuffer(obj)
% Reinitialize observation buffer. Allocate as dlarray to
% support automatic differentiation with dlfeval and
% dlgradient.

%format：CBT

obj.obsBuffer = dlarray(…
zeros(obj.numObs,obj.options.MaxBufferLen),’CB’);

% Reinitialize action buffer with valid actions.
obj.actionBuffer = dlarray(…
zeros(obj.numAct,obj.options.MaxBufferLen),’CB’);

% Reinitialize reward buffer.
obj.rewardBuffer = dlarray(zeros(1,obj.options.MaxBufferLen),’CB’);

% Reinitialize nextState buffer.
obj.nextObsBuffer = dlarray(…
zeros(obj.numObs,obj.options.MaxBufferLen),’CB’);

% Reinitialize mask buffer.
obj.isDoneBuffer = dlarray(zeros(1,obj.options.MaxBufferLen),’CB’);

end

%Create networks
%Actor
function actor = CreateActor(obj,numObs,numAct,hid_dim,obsInfo,actInfo)
% Create the actor network layers.
commonPath = [
featureInputLayer(numObs,Name="obsInLyr")
fullyConnectedLayer(hid_dim)
layerNormalizationLayer
reluLayer
fullyConnectedLayer(hid_dim)
layerNormalizationLayer
reluLayer(Name="comPathOutLyr")
];
meanPath = [
fullyConnectedLayer(numAct,Name="meanOutLyr")
];
stdPath = [
fullyConnectedLayer(numAct,Name="stdInLyr")
softplusLayer(Name="stdOutLyr")
];
% Connect the layers.
actorNetwork = layerGraph(commonPath);
actorNetwork = addLayers(actorNetwork,meanPath);
actorNetwork = addLayers(actorNetwork,stdPath);
actorNetwork = connectLayers(actorNetwork,"comPathOutLyr","meanOutLyr/in");
actorNetwork = connectLayers(actorNetwork,"comPathOutLyr","stdInLyr/in");

actordlnet = dlnetwork(actorNetwork);
actor = initialize(actordlnet);
end

%Critic
function [critic1,critic2,critic_target1,critic_target2] = CreateCritic(obj,numObs,numAct,hid_dim,obsInfo,actInfo)
% Define the network layers.
criticNet = [
featureInputLayer(numObs+numAct,Name="obsInLyr")%input:[obs act]
fullyConnectedLayer(hid_dim)
layerNormalizationLayer
reluLayer
fullyConnectedLayer(hid_dim)
layerNormalizationLayer
reluLayer
fullyConnectedLayer(1,Name="QValueOutLyr")
];

% Connect the layers.
criticNet = layerGraph(criticNet);
criticDLnet = dlnetwork(criticNet,’Initialize’,false);
critic1 = initialize(criticDLnet);
critic2 = initialize(criticDLnet);%c1 and c2 different initilization
critic_target1 = initialize(criticDLnet);
critic_target1.Learnables = critic1.Learnables;
critic_target1.State = critic1.State;
critic_target2 = initialize(criticDLnet);
critic_target2.Learnables = critic2.Learnables;
critic_target2.State = critic2.State;

end

function logP = logProbBoundedAction(obj,boundedAction,mu,sigma)
%used to calculate log probability for tanh(gaussian)
%validated, nothing wrong with this function
eps=1e-10;
logP = sum(log(1/sqrt(2*pi)./sigma.*exp(-0.5*(0.5*…
log((1+boundedAction+eps)./(1-boundedAction+eps))-mu).^2./sigma.^2).*1./(1-boundedAction.^2+eps)),1);
end

%loss functions
function [vLoss_1, vLoss_2, criticGrad_1, criticGrad_2] = criticLoss(obj,batchExperiences,c1,c2)
batchObs = batchExperiences{1};
batchAction = batchExperiences{2};
batchReward = batchExperiences{3};
batchNextObs = batchExperiences{4};
batchIsDone = batchExperiences{5};
batchSize = size(batchObs,2);

gamma = obj.options.DiscountFactor;
y = dlarray(zeros(1,batchSize));%CB(C=1)
y = y + batchReward;
actionNext = getActionWithExploration_dlarray(obj,batchNextObs);%CB
actionNext = actionNext{1};
Qt1=predict(obj.critic_target1,cat(1,batchNextObs,actionNext));%CB(C=1)
Qt2=predict(obj.critic_target2,cat(1,batchNextObs,actionNext));%CB(C=1)

[mu,sigma] = predict(obj.actor,batchNextObs);%CB:numAct*batch
next_action = tanh(mu + sigma.*randn(size(sigma)));
logP = logProbBoundedAction(obj,next_action,mu,sigma);
y = y + (1 – batchIsDone).*(gamma*(min(cat(1,Qt1,Qt2),[],1) – exp(obj.log_alpha)*logP));

critic_input = cat(1,batchObs,batchAction);
Q1 = forward(c1,critic_input);
Q2 = forward(c2,critic_input);
vLoss_1 = 1/2*mean((y – Q1).^2,’all’);
vLoss_2 = 1/2*mean((y – Q2).^2,’all’);
criticGrad_1 = dlgradient(vLoss_1,c1.Learnables);
criticGrad_2 = dlgradient(vLoss_2,c2.Learnables);
end

function [aLoss,actorGrad] = actorLoss(obj,batchExperiences,actor)
batchObs = batchExperiences{1};
batchSize = size(batchObs,2);
[mu,sigma] = forward(actor,batchObs);%CB:numAct*batch
curr_action = tanh(mu + sigma.*randn(size(sigma)));%reparameterization
critic_input = cat(1,batchObs,curr_action);
Q1=forward(obj.critic1,critic_input);%CB(C=1)
Q2=forward(obj.critic2,critic_input);%CB(C=1)

logP = logProbBoundedAction(obj,curr_action,mu,sigma);
aLoss = mean(-min(cat(1,Q1,Q2),[],1) + exp(obj.log_alpha) * logP,’all’);

actorGrad= dlgradient(aLoss,actor.Learnables);
end

function [eLoss,entGrad] = entropyLoss(obj,batchExperiences,logAlpha)
batchObs = batchExperiences{1};
[mu,sigma] = predict(obj.actor,batchObs);%CB:numAct*batch
curr_action = tanh(mu + sigma.*randn(size(sigma)));
ent = mean(-logProbBoundedAction(obj,curr_action,mu,sigma));
eLoss = exp(logAlpha) * (ent – obj.options.TargetEntropy);
entGrad = dlgradient(eLoss,logAlpha);
end

end
methods(Access=protected)
%return SampleTime
function ts = getSampleTime_(obj)
ts = obj.Ts;
end
%get action without exploration
function action = getActionImpl(obj,obs)
%obs:dlarray CB
if ~isa(obs,’dlarray’)
if isa(obs,’cell’)
obs = dlarray(obs{1},’CB’);
else
obs = dlarray(obs,’CB’);
end
end
[mu,~] = predict(obj.actor,obs);
mu = extractdata(mu);
action = {tanh(mu)};
end

%get action with exploration
function action = getActionWithExplorationImpl(obj,obs)
%obs:dlarray CT
if ~isa(obs,’dlarray’) || size(obs,1)~=obj.numObs
obs = dlarray(randn(obj.numObs,1),’CB’);
end
[mu,sigma] = predict(obj.actor,obs);
mu = extractdata(mu);
sigma = extractdata(sigma);
action = {tanh(mu + sigma .* randn(size(sigma)))};
end

function action = getActionWithExploration_dlarray(obj,obs)
[mu,sigma] = predict(obj.actor,obs);
action = {tanh(mu + sigma .* randn(size(sigma)))};
end

%learning
function action = learnImpl(obj,Experience)
% Extract data from experience.

obs = Experience{1};
action = Experience{2};
reward = Experience{3};
nextObs = Experience{4};
isDone = logical(Experience{5});

obj.obsBuffer(:,obj.bufferIdx+1,:) = obs{1};
obj.actionBuffer(:,obj.bufferIdx+1,:) = action{1};
obj.rewardBuffer(:,obj.bufferIdx+1) = reward;
obj.nextObsBuffer(:,obj.bufferIdx+1,:) = nextObs{1};
obj.isDoneBuffer(:,obj.bufferIdx+1) = isDone;
obj.bufferLen = max(obj.bufferLen,obj.bufferIdx+1);
obj.bufferIdx = mod(obj.bufferIdx+1,obj.options.MaxBufferLen);

if obj.bufferLen>=max(obj.options.WarmUpSteps,obj.options.MiniBatchSize)
obj.counter = obj.counter + 1;
if (obj.options.LearningFrequency==-1 && isDone) || …
(obj.options.LearningFrequency>0 && mod(obj.counter,obj.options.LearningFrequency)==0)
for gstep = 1:obj.options.NumGradientStepsPerUpdate
%sample batch
batchSize = obj.options.MiniBatchSize;
batchInd = randperm(obj.bufferLen,batchSize);

batchExperience = {
obj.obsBuffer(:,batchInd,:),…
obj.actionBuffer(:,batchInd,:),…
obj.rewardBuffer(:,batchInd),…
obj.nextObsBuffer(:,batchInd,:),…
obj.isDoneBuffer(:,batchInd)
};

%update the parameters of each critic
[cLoss1,cLoss2,criticGrad_1,criticGrad_2] = dlfeval(@(x,c1,c2)obj.criticLoss(x,c1,c2),batchExperience,obj.critic1,obj.critic2);
obj.cLoss = min(extractdata(cLoss1),extractdata(cLoss2));
[obj.critic1.Learnables.Value,obj.criticOptimizer_1] = update(obj.criticOptimizer_1,obj.critic1.Learnables.Value,criticGrad_1.Value);
[obj.critic2.Learnables.Value,obj.criticOptimizer_2] = update(obj.criticOptimizer_2,obj.critic2.Learnables.Value,criticGrad_2.Value);

if (mod(obj.counter,obj.options.PolicyUpdateFrequency)==0 && obj.options.LearningFrequency==-1) ||…
(mod(obj.counter,obj.options.LearningFrequency * obj.options.PolicyUpdateFrequency)==0 …
&& obj.options.LearningFrequency>0)
%update the parameters of actor
[aloss,actorGrad] = dlfeval(…
@(x,actor)obj.actorLoss(x,actor),…
batchExperience,obj.actor);
obj.aLoss = extractdata(aloss);
[obj.actor.Learnables.Value,obj.actorOptimizer] = update(obj.actorOptimizer,obj.actor.Learnables.Value,actorGrad.Value);

%update the entropy weight
[eloss,entGrad] = dlfeval(@(x,alpha)obj.entropyLoss(x,alpha),batchExperience,obj.log_alpha);
obj.eLoss = extractdata(eloss);
% disp(obj.alpha)
[obj.log_alpha,obj.entWgtOptimizer] = update(obj.entWgtOptimizer,{obj.log_alpha},{entGrad});
obj.log_alpha = obj.log_alpha{1};
end

%update critic targets
%1
critic1_params = obj.critic1.Learnables.Value;%cell array network params
critic_target1_params = obj.critic_target1.Learnables.Value;
for i=1:size(critic1_params,1)
obj.critic_target1.Learnables.Value{i} = obj.options.TargetSmoothFactor * critic1_params{i}…
+ (1 – obj.options.TargetSmoothFactor) * critic_target1_params{i};
end
%2
critic2_params = obj.critic2.Learnables.Value;%cell array network params
critic_target2_params = obj.critic_target2.Learnables.Value;
for i=1:size(critic2_params,1)
obj.critic_target2.Learnables.Value{i} = obj.options.TargetSmoothFactor * critic2_params{i}…
+ (1 – obj.options.TargetSmoothFactor) * critic_target2_params{i};
end
% end
end
end

end
action = getActionWithExplorationImpl(obj,nextObs{1});
end
end
end

2.Configuration of ‘options’ property(same as those used for the built-in SAC agent)
options.MaxBufferLen = 1e4;
options.WarmUpSteps = 1000;
options.MiniBatchSize = 256;
options.LearningFrequency = -1;%when -1: train after each episode
options.EntropyLossWeight = 1;
options.DiscountFactor = 0.99;
options.PolicyUpdateFrequency = 1;
options.TargetEntropy = -2;
options.TargetUpdateFrequency = 1;
options.TargetSmoothFactor = 1e-3;
options.NumGradientStepsPerUpdate = 10;
%optimizerOptions: actor critic1 critic2 entWgt(alpha)
%encoder decoder
options.OptimizerOptions = {
rlOptimizerOptions("Algorithm","adam","GradientThreshold",1,’LearnRate’,1e-3),…
rlOptimizerOptions("Algorithm","adam","GradientThreshold",1,’LearnRate’,1e-3),…
rlOptimizerOptions("Algorithm","adam","GradientThreshold",1,’LearnRate’,1e-3),…
rlOptimizerOptions("Algorithm","adam",’LearnRate’,3e-4),…
rlOptimizerOptions("Algorithm","adam","GradientThreshold",1,’LearnRate’,1e-3),…
rlOptimizerOptions("Algorithm","adam","GradientThreshold",1,’LearnRate’,1e-3)};
options.base_seed=940;
3.training

clc;
clear;
close all;
run(‘init_car_params.m’);
%create RL env
numObs = 4; % vx vy r beta_user
numAct = 2; % st_angle_ref rw_omega_ref

obsInfo = rlNumericSpec([numObs 1]);

actInfo = rlNumericSpec([numAct 1]);
actInfo.LowerLimit = -1;
actInfo.UpperLimit = 1;

mdl = "prius_sm_model";
blk = mdl + "/RL Agent";
env = rlSimulinkEnv(mdl,blk,obsInfo,actInfo);

params=struct(‘rw_radius’,rw_radius,’a’,a,’b’,b,’init_vx’,init_vx,’init_yaw_rate’,init_yaw_rate);
env.ResetFcn = @(in) PriusResetFcn(in,params,mdl);

Ts = 1/10;
Tf = 5;

%create actor
rnd_seed=940;
algorithm = ‘MySAC’;
switch algorithm
case ‘SAC’
agent = createNetworks(rnd_seed,numObs,numAct,obsInfo,actInfo,Ts);
case ‘MySAC’
hid_dim = 256;
options=getDWMLAgentOptions();
agent = MySACAgent(numObs,numAct,obsInfo,actInfo,hid_dim,Ts,options);
end

%%
%train agent
close all
maxEpisodes = 6000;
maxSteps = floor(Tf/Ts);
useParallel = false;
run_idx=9;
saveAgentDir = [‘savedAgents/’,algorithm,’/’,num2str(run_idx)];
switch algorithm
case ‘SAC’
trainOpts = rlTrainingOptions(…
MaxEpisodes=maxEpisodes, …
MaxStepsPerEpisode=maxSteps, …
ScoreAveragingWindowLength=100, …
Plots="training-progress", …
StopTrainingCriteria="AverageReward", …
UseParallel=useParallel,…
SaveAgentCriteria=’EpisodeReward’,…
SaveAgentValue=35,…
SaveAgentDirectory=saveAgentDir);
% SaveAgentCriteria=’EpisodeFrequency’,…
% SaveAgentValue=1,…
case ‘MySAC’
trainOpts = rlTrainingOptions(…
MaxEpisodes=maxEpisodes, …
MaxStepsPerEpisode=maxSteps, …
ScoreAveragingWindowLength=100, …
Plots="training-progress", …
StopTrainingCriteria="AverageReward", …
UseParallel=useParallel,…
SaveAgentCriteria=’EpisodeReward’,…
SaveAgentValue=35,…
SaveAgentDirectory=saveAgentDir);
end

set_param(mdl,"FastRestart","off");%for random initialization
if trainOpts.UseParallel
% Disable visualization in Simscape Mechanics Explorer
set_param(mdl, SimMechanicsOpenEditorOnUpdate="off");
save_system(mdl);
else
% Enable visualization in Simscape Mechanics Explorer
set_param(mdl, SimMechanicsOpenEditorOnUpdate="on");
end

%load training data
monitor = trainingProgressMonitor();
logger = rlDataLogger(monitor);
logger.EpisodeFinishedFcn = @myEpisodeLoggingFcn;

doTraining = true;
if doTraining
trainResult = train(agent,env,trainOpts,Logger=logger);
end

% %logger callback used for MySACAgent
function dataToLog = myEpisodeLoggingFcn(data)
dataToLog.criticLoss = data.Agent.cLoss;
dataToLog.actorLoss = data.Agent.aLoss;
dataToLog.entLoss = data.Agent.eLoss;
% dataToLog.denoiseLoss = data.Agent.dnLoss;
end
In the simulink environment used, action output by the Agent block(in [-1,1]) is denormalized and fed into the environment.

I think possible causes of the problem include:
1.Wrong implementation of critic loss. As shown in the training progress, critic loss seemed to diverge. It’s hardly caused by hyperparameters(batch size or learning rate or target update frequency) because they worked well for the built-in agent. So it is more likely the critic loss is wrong.
2.Wrong implementation of replay buffer. I implemented the replay buffer as a circular queue, where I sampled uniformly to get batch training data. From the comparison of the training progress shown above, the custom SAC agent did explore states with high reward(around 30) but failed to exploit them, So I guess there is still problem with my replay buffer.
3.Gradient flow was broken.The learning is done with the help of MATLAB deep learning automatic differentiation. Perhaps some of my implementation violates the computational rule of automatic differentiation, which broke the gradient flow during forward computation or backpropagation and led to wrong result.
4.Gradient step(update frequency). In current implementation, NumGradientStepsPerUpdate gradient steps are executed after each episode. During each gradient step, cirtic(s) and actor, alongwith entropy weight, is updated once. I am not sure whether the current implementation of gradient step has got the update frequency right.
5.Also could be normalization problem, but I am not so sure.
I plan to debug 3 first.
Please read the code and help find potential causes of the gap between the custom SAC agent and the built-in one.
Finally, I am actually trying to extend SAC algorithm to a more complex framework. I didn’t choose to inherit the built-in SAC agent(rlSACAgent), would it be recommended to do my development by doing so? reinforcement learning, soft actor critic, sac, automatic differentiation, custom agent MATLAB Answers — New Questions

​

The following is part of the program. I want to know how to extract the weight values from the trained DQN network.

DQNnet = [
imageInputLayer([1 520 1],"Name","ImageFeatureInput","Normalization","none")
fullyConnectedLayer(1024,"Name","fc1")
reluLayer("Name","relu1")
% fullyConnectedLayer(512,"Name","fc2")
% reluLayer("Name","relu2")
fullyConnectedLayer(14,"Name","fc3")
softmaxLayer("Name","softmax")
classificationLayer("Name","ActionOutput")];

ObsInfo = getObservationInfo(env);
ActInfo = getActionInfo(env);

DQNOpts = rlRepresentationOptions(‘LearnRate’,0.0001,’GradientThreshold’,1,’UseDevice’,’gpu’);
DQNagent = rlQValueRepresentation(DQNnet,ObsInfo,ActInfo,’Observation’,{‘ImageFeatureInput’},’ActionInputNames’,{‘BoundingBox Actions’},DQNOpts);

agentOpts = rlDQNAgentOptions(…
‘UseDoubleDQN’,true …
,’MiniBatchSize’,256);

agentOpts.EpsilonGreedyExploration.Epsilon = 1;

agent = rlDQNAgent(DQNagent,agentOpts);

%% Agent Training
% Training options

trainOpts = rlTrainingOptions(…
‘MaxEpisodes’, 100, …
‘MaxStepsPerEpisode’, 100, …
‘Verbose’, true, …
‘Plots’,’training-progress’,…
‘ScoreAveragingWindowLength’,400,…
‘StopTrainingCriteria’,’AverageSteps’,…
‘StopTrainingValue’,1000000000,…
‘SaveAgentDirectory’, pwd + "agents");

% Agent training
trainingStats = train(agent,env,trainOpts);The following is part of the program. I want to know how to extract the weight values from the trained DQN network.

DQNnet = [
imageInputLayer([1 520 1],"Name","ImageFeatureInput","Normalization","none")
fullyConnectedLayer(1024,"Name","fc1")
reluLayer("Name","relu1")
% fullyConnectedLayer(512,"Name","fc2")
% reluLayer("Name","relu2")
fullyConnectedLayer(14,"Name","fc3")
softmaxLayer("Name","softmax")
classificationLayer("Name","ActionOutput")];

ObsInfo = getObservationInfo(env);
ActInfo = getActionInfo(env);

DQNOpts = rlRepresentationOptions(‘LearnRate’,0.0001,’GradientThreshold’,1,’UseDevice’,’gpu’);
DQNagent = rlQValueRepresentation(DQNnet,ObsInfo,ActInfo,’Observation’,{‘ImageFeatureInput’},’ActionInputNames’,{‘BoundingBox Actions’},DQNOpts);

agentOpts = rlDQNAgentOptions(…
‘UseDoubleDQN’,true …
,’MiniBatchSize’,256);

agentOpts.EpsilonGreedyExploration.Epsilon = 1;

agent = rlDQNAgent(DQNagent,agentOpts);

%% Agent Training
% Training options

trainOpts = rlTrainingOptions(…
‘MaxEpisodes’, 100, …
‘MaxStepsPerEpisode’, 100, …
‘Verbose’, true, …
‘Plots’,’training-progress’,…
‘ScoreAveragingWindowLength’,400,…
‘StopTrainingCriteria’,’AverageSteps’,…
‘StopTrainingValue’,1000000000,…
‘SaveAgentDirectory’, pwd + "agents");

% Agent training
trainingStats = train(agent,env,trainOpts); The following is part of the program. I want to know how to extract the weight values from the trained DQN network.

DQNnet = [
imageInputLayer([1 520 1],"Name","ImageFeatureInput","Normalization","none")
fullyConnectedLayer(1024,"Name","fc1")
reluLayer("Name","relu1")
% fullyConnectedLayer(512,"Name","fc2")
% reluLayer("Name","relu2")
fullyConnectedLayer(14,"Name","fc3")
softmaxLayer("Name","softmax")
classificationLayer("Name","ActionOutput")];

ObsInfo = getObservationInfo(env);
ActInfo = getActionInfo(env);

DQNOpts = rlRepresentationOptions(‘LearnRate’,0.0001,’GradientThreshold’,1,’UseDevice’,’gpu’);
DQNagent = rlQValueRepresentation(DQNnet,ObsInfo,ActInfo,’Observation’,{‘ImageFeatureInput’},’ActionInputNames’,{‘BoundingBox Actions’},DQNOpts);

agentOpts = rlDQNAgentOptions(…
‘UseDoubleDQN’,true …
,’MiniBatchSize’,256);

agentOpts.EpsilonGreedyExploration.Epsilon = 1;

agent = rlDQNAgent(DQNagent,agentOpts);

%% Agent Training
% Training options

trainOpts = rlTrainingOptions(…
‘MaxEpisodes’, 100, …
‘MaxStepsPerEpisode’, 100, …
‘Verbose’, true, …
‘Plots’,’training-progress’,…
‘ScoreAveragingWindowLength’,400,…
‘StopTrainingCriteria’,’AverageSteps’,…
‘StopTrainingValue’,1000000000,…
‘SaveAgentDirectory’, pwd + "agents");

% Agent training
trainingStats = train(agent,env,trainOpts); drl, neural network, dqn MATLAB Answers — New Questions

​

The code I have is where f is a function handle, a is a real number, and n is a positive integer:

function r=mynewton(f,a,n)
syms x
f=@x;
c=f(x);
y(1)=a;
for i=[1:length(n)]
y(i+1)=y(i)-(c(i)/diff(c(i)));
end;
r=y

The erros I am getting are:

Undefined function or variable ‘x’.
Error in mynewton (line 4)
c=f(x);

How do I fix these errors?The code I have is where f is a function handle, a is a real number, and n is a positive integer:

function r=mynewton(f,a,n)
syms x
f=@x;
c=f(x);
y(1)=a;
for i=[1:length(n)]
y(i+1)=y(i)-(c(i)/diff(c(i)));
end;
r=y

The erros I am getting are:

Undefined function or variable ‘x’.
Error in mynewton (line 4)
c=f(x);

How do I fix these errors? The code I have is where f is a function handle, a is a real number, and n is a positive integer:

function r=mynewton(f,a,n)
syms x
f=@x;
c=f(x);
y(1)=a;
for i=[1:length(n)]
y(i+1)=y(i)-(c(i)/diff(c(i)));
end;
r=y

The erros I am getting are:

Undefined function or variable ‘x’.
Error in mynewton (line 4)
c=f(x);

How do I fix these errors? newton raphson method function handle MATLAB Answers — New Questions

​

Hi everyone! I am working on a project where I need to calibrate multiple cameras observing a scene, to ultimately be able to get 3D points of an object in later videos collected by the same cameras. The cameras are stationary. Importantly, I need to be able to triangulate the checkerboard points from the calibration and then do sparse bundle adjustment on these points to improve the acuracy of the camera pose estimation and 3D checkerboard points from the calibration. Sparse bundle adjustment (bundleAdjustment) can take in either pointTrack objects or worldpointset objects.
I have two calibration sessions (front camera and rear right, and the front camera and rear left – they are in a triangular config) from which I load the stereoParams, I have also stored the useful data in a structure called ‘s’.
I then get the 3D coordinates of the checkerboards, and using the ‘worldpointset’ and feature matching approach. I have included all code (including the code I used to save important variables).
The error I get with the bundleAdjustment function is the following:
Error using vision.internal.bundleAdjust.validateAndParseInputs
The number of feature points in view 1 must be greater than or equal to 51.
Error in vision.internal.bundleAdjust.sparseBA (line 39)
vision.internal.bundleAdjust.validateAndParseInputs(optimType, mfilename, varargin{:});
Error in bundleAdjustment (line 10)
vision.internal.bundleAdjust.sparseBA(‘full’, mfilename, varargin{:});

When I investigated using pointTrack, it seems that it is best used for tracking a point through multiple frames in a video, but not great for my application where I want to track one point through 3 different views at once.

AT LAST–> MY QUESTION:
Am I using worldpointset correctly for this application, and if so, can someone please help me figure out where this error in the feature points is coming from?
If not, would pointTrack be better for my application if I change the dimensionality of my problem? If pointTrack is better, I would need to track a point through the frames of each camera and somehow correlate and triangulate points that way.
**Note, with structure ‘s’ containing many images it was too large to upload (even when compressed), so I uploaded a screenshot of the structure. But hopefully my code helps with context. The visualisation runs though!
load("params.mat","params")
intrinsics1 = params.cam1.Intrinsics;
intrinsics2 = params.cam2.Intrinsics;
intrinsics3 = params.cam3.Intrinsics;
intrinsics4 = params.cam4.Intrinsics;

intrinsicsFront = intrinsics2;
intrinsicsRLeft = intrinsics3;
intrinsicsRRight = intrinsics4;

%% Visualise cameras
load("stereoParams1.mat")
load("stereoParams2.mat")

figure; showExtrinsics(stereoParams1, ‘CameraCentric’)
hold on;
showExtrinsics(stereoParams2, ‘CameraCentric’);
hold off;

%initialise camera 1 pose as at 0, with no rotation
front_absolute_pose = rigidtform3d([0 0 0], [0 0 0]);

%% Get 3D Points
load("s_struct.mat","s")

board_shape = s.front.board_shape;
camPoseVSet = imageviewset;
camPoseVSet = addView(camPoseVSet,1,front_absolute_pose);
camPoseVSet = addView(camPoseVSet,2,stereoParams1.PoseCamera2);
camPoseVSet = addView(camPoseVSet,3,stereoParams2.PoseCamera2);
camposes = poses(camPoseVSet);
intrinsicsArray = [intrinsicsFront, intrinsicsRRight, intrinsicsRLeft];
frames = fieldnames(s.front.points);
framesRearRight = fieldnames(s.rearright.points);
framesRearLeft = fieldnames(s.rearleft.points);
wpSet =worldpointset;
wpsettrial = worldpointset;
for i =1:length(frames) %for frames in front
frame_i = frames{i};
pointsFront = s.front.points.(frame_i);
pointsFrontUS = s.front.unshapedPoints.(frame_i);
container = contains(framesRearRight, frame_i);

j = 1;
if ~isempty(container) && any(container)
pointsRearRight = s.rearright.points.(frame_i);
pointsRearRightUS = s.rearright.unshapedPoints.(frame_i);
pointIn3D = [];
pointIn3Dnew = [];
[features1, validPts1] = extractFeatures(im2gray(s.front.imageFile.(frame_i)), pointsFrontUS);
[features2, validPts2] = extractFeatures(im2gray(s.rearright.imageFile.(frame_i)), pointsRearRightUS);
indexPairs = matchFeatures(features1,features2);
matchedPoints1 = validPts1(indexPairs(:,1),:);
matchedPoints2 = validPts2(indexPairs(:,2),:);
worldPTS = triangulate(matchedPoints1, matchedPoints2, stereoParams2);
[wpsettrial,newPointIndices] = addWorldPoints(wpsettrial,worldPTS);
wpsettrial = addCorrespondences(wpsettrial,1,newPointIndices,indexPairs(:,1));
wpsettrial = addCorrespondences(wpsettrial,3,newPointIndices,indexPairs(:,2));
sz = size(s.front.points.(frame_i));

for h =1: sz(1)

for w = 1:sz(2)
point2track = [pointsFront(h,w,1), pointsFront(h,w,2); pointsRearRight(h,w,1), pointsRearRight(h,w,2)];
IDs = [1, 3];
track = pointTrack(IDs,point2track);
triang3D = triangulate([pointsFront(h,w,1), pointsFront(h,w,2)], [pointsRearRight(h,w,1), pointsRearRight(h,w,2)], stereoParams1);
% [wpSet,newPointIndices] = addWorldPoints(wpSet,triang3D);
% wpSet = addCorrespondences(wpSet,1,j,j);
% wpSet = addCorrespondences(wpSet,3,j,j);
pointIn3D = [pointIn3D;triang3D];
j=j+1;

end

end

pointIn3D = reshape3D(pointIn3D, board_shape);
%xyzPoints = reshape3D(pointIn3D,board_shape);
s.frontANDright.PT3D.(frame_i) = pointIn3D;
%s.frontANDright.PT3DSBA.(frame_i) = xyzPoints;
end
container = contains(framesRearLeft, frame_i);
m=1;
if ~isempty(container) && any(container)

pointsRearLeft = s.rearleft.points.(frame_i);
pointsRearLeftUS = s.rearleft.unshapedPoints.(frame_i);
pointIn3D = [];
pointIn3Dnew = [];
sz = size(s.front.points.(frame_i));
[features1, validPts1] = extractFeatures(im2gray(s.front.imageFile.(frame_i)), pointsFrontUS);
[features2, validPts2] = extractFeatures(im2gray(s.rearleft.imageFile.(frame_i)), pointsRearLeftUS);
indexPairs = matchFeatures(features1,features2);
matchedPoints1 = validPts1(indexPairs(:,1),:);
matchedPoints2 = validPts2(indexPairs(:,2),:);
worldPTS = triangulate(matchedPoints1, matchedPoints2, stereoParams1);
[wpsettrial,newPointIndices] = addWorldPoints(wpsettrial,worldPTS);
wpsettrial = addCorrespondences(wpsettrial,1,newPointIndices,indexPairs(:,1));
wpsettrial = addCorrespondences(wpsettrial,2,newPointIndices,indexPairs(:,2));

for h =1: sz(1)

for w = 1:sz(2)
point2track = [pointsFront(h,w,1), pointsFront(h,w,2); pointsRearLeft(h,w,1), pointsRearLeft(h,w,2)];

IDs = [1, 2];
track = pointTrack(IDs,point2track);
triang3D = triangulate([pointsFront(h,w,1), pointsFront(h,w,2)], [pointsRearLeft(h,w,1), pointsRearLeft(h,w,2)], stereoParams1);
% wpSet = addWorldPoints(wpSet,triang3D);
% wpSet = addCorrespondences(wpSet,1,m,m);
% wpSet = addCorrespondences(wpSet,2,m,m);
pointIn3D = [pointIn3D;triang3D];
m = m+1;

end

end
pointIn3D = reshape3D(pointIn3D, board_shape);
%xyzPoints = reshape3D(pointIn3D,board_shape);
s.frontANDleft.PT3D.(frame_i) = pointIn3D;
%s.frontANDleft.PT3DSBA.(frame_i) = xyzPoints;
end
[wpSetRefined,vSetRefined,pointIndex] = bundleAdjustment(wpsettrial,camPoseVSet,[1,3,2],intrinsicsArray, FixedViewIDs=[1,3,2], …
Solver="preconditioned-conjugate-gradient")
end

function [img_name, ptsUS,pts, worldpoints] = reformatData(img_name, pts, board_shape, worldpoints)
%method taken from acinoset code
img_name = img_name(1:strfind(img_name,’ ‘)-1);
img_name = replace(img_name, ‘.’,’_’);
ptsUS = pts;
pts = pagetranspose(reshape(pts, [board_shape, 2]));
pts = pagetranspose(reshape(pts, [board_shape, 2])); %repetition is purposeful

worldpoints = pagetranspose(reshape(worldpoints, [board_shape,2]));
worldpoints = pagetranspose(reshape(worldpoints, [board_shape,2]));
end

function pts = reshape3D(points3D, board_shape)
pts = pagetranspose(reshape(points3D, [board_shape, 3]));
pts = pagetranspose(reshape(pts, [board_shape, 3])); %repetition is purposeful

endHi everyone! I am working on a project where I need to calibrate multiple cameras observing a scene, to ultimately be able to get 3D points of an object in later videos collected by the same cameras. The cameras are stationary. Importantly, I need to be able to triangulate the checkerboard points from the calibration and then do sparse bundle adjustment on these points to improve the acuracy of the camera pose estimation and 3D checkerboard points from the calibration. Sparse bundle adjustment (bundleAdjustment) can take in either pointTrack objects or worldpointset objects.
I have two calibration sessions (front camera and rear right, and the front camera and rear left – they are in a triangular config) from which I load the stereoParams, I have also stored the useful data in a structure called ‘s’.
I then get the 3D coordinates of the checkerboards, and using the ‘worldpointset’ and feature matching approach. I have included all code (including the code I used to save important variables).
The error I get with the bundleAdjustment function is the following:
Error using vision.internal.bundleAdjust.validateAndParseInputs
The number of feature points in view 1 must be greater than or equal to 51.
Error in vision.internal.bundleAdjust.sparseBA (line 39)
vision.internal.bundleAdjust.validateAndParseInputs(optimType, mfilename, varargin{:});
Error in bundleAdjustment (line 10)
vision.internal.bundleAdjust.sparseBA(‘full’, mfilename, varargin{:});

When I investigated using pointTrack, it seems that it is best used for tracking a point through multiple frames in a video, but not great for my application where I want to track one point through 3 different views at once.

AT LAST–> MY QUESTION:
Am I using worldpointset correctly for this application, and if so, can someone please help me figure out where this error in the feature points is coming from?
If not, would pointTrack be better for my application if I change the dimensionality of my problem? If pointTrack is better, I would need to track a point through the frames of each camera and somehow correlate and triangulate points that way.
**Note, with structure ‘s’ containing many images it was too large to upload (even when compressed), so I uploaded a screenshot of the structure. But hopefully my code helps with context. The visualisation runs though!
load("params.mat","params")
intrinsics1 = params.cam1.Intrinsics;
intrinsics2 = params.cam2.Intrinsics;
intrinsics3 = params.cam3.Intrinsics;
intrinsics4 = params.cam4.Intrinsics;

intrinsicsFront = intrinsics2;
intrinsicsRLeft = intrinsics3;
intrinsicsRRight = intrinsics4;

%% Visualise cameras
load("stereoParams1.mat")
load("stereoParams2.mat")

figure; showExtrinsics(stereoParams1, ‘CameraCentric’)
hold on;
showExtrinsics(stereoParams2, ‘CameraCentric’);
hold off;

%initialise camera 1 pose as at 0, with no rotation
front_absolute_pose = rigidtform3d([0 0 0], [0 0 0]);

%% Get 3D Points
load("s_struct.mat","s")

board_shape = s.front.board_shape;
camPoseVSet = imageviewset;
camPoseVSet = addView(camPoseVSet,1,front_absolute_pose);
camPoseVSet = addView(camPoseVSet,2,stereoParams1.PoseCamera2);
camPoseVSet = addView(camPoseVSet,3,stereoParams2.PoseCamera2);
camposes = poses(camPoseVSet);
intrinsicsArray = [intrinsicsFront, intrinsicsRRight, intrinsicsRLeft];
frames = fieldnames(s.front.points);
framesRearRight = fieldnames(s.rearright.points);
framesRearLeft = fieldnames(s.rearleft.points);
wpSet =worldpointset;
wpsettrial = worldpointset;
for i =1:length(frames) %for frames in front
frame_i = frames{i};
pointsFront = s.front.points.(frame_i);
pointsFrontUS = s.front.unshapedPoints.(frame_i);
container = contains(framesRearRight, frame_i);

j = 1;
if ~isempty(container) && any(container)
pointsRearRight = s.rearright.points.(frame_i);
pointsRearRightUS = s.rearright.unshapedPoints.(frame_i);
pointIn3D = [];
pointIn3Dnew = [];
[features1, validPts1] = extractFeatures(im2gray(s.front.imageFile.(frame_i)), pointsFrontUS);
[features2, validPts2] = extractFeatures(im2gray(s.rearright.imageFile.(frame_i)), pointsRearRightUS);
indexPairs = matchFeatures(features1,features2);
matchedPoints1 = validPts1(indexPairs(:,1),:);
matchedPoints2 = validPts2(indexPairs(:,2),:);
worldPTS = triangulate(matchedPoints1, matchedPoints2, stereoParams2);
[wpsettrial,newPointIndices] = addWorldPoints(wpsettrial,worldPTS);
wpsettrial = addCorrespondences(wpsettrial,1,newPointIndices,indexPairs(:,1));
wpsettrial = addCorrespondences(wpsettrial,3,newPointIndices,indexPairs(:,2));
sz = size(s.front.points.(frame_i));

for h =1: sz(1)

for w = 1:sz(2)
point2track = [pointsFront(h,w,1), pointsFront(h,w,2); pointsRearRight(h,w,1), pointsRearRight(h,w,2)];
IDs = [1, 3];
track = pointTrack(IDs,point2track);
triang3D = triangulate([pointsFront(h,w,1), pointsFront(h,w,2)], [pointsRearRight(h,w,1), pointsRearRight(h,w,2)], stereoParams1);
% [wpSet,newPointIndices] = addWorldPoints(wpSet,triang3D);
% wpSet = addCorrespondences(wpSet,1,j,j);
% wpSet = addCorrespondences(wpSet,3,j,j);
pointIn3D = [pointIn3D;triang3D];
j=j+1;

end

end

pointIn3D = reshape3D(pointIn3D, board_shape);
%xyzPoints = reshape3D(pointIn3D,board_shape);
s.frontANDright.PT3D.(frame_i) = pointIn3D;
%s.frontANDright.PT3DSBA.(frame_i) = xyzPoints;
end
container = contains(framesRearLeft, frame_i);
m=1;
if ~isempty(container) && any(container)

pointsRearLeft = s.rearleft.points.(frame_i);
pointsRearLeftUS = s.rearleft.unshapedPoints.(frame_i);
pointIn3D = [];
pointIn3Dnew = [];
sz = size(s.front.points.(frame_i));
[features1, validPts1] = extractFeatures(im2gray(s.front.imageFile.(frame_i)), pointsFrontUS);
[features2, validPts2] = extractFeatures(im2gray(s.rearleft.imageFile.(frame_i)), pointsRearLeftUS);
indexPairs = matchFeatures(features1,features2);
matchedPoints1 = validPts1(indexPairs(:,1),:);
matchedPoints2 = validPts2(indexPairs(:,2),:);
worldPTS = triangulate(matchedPoints1, matchedPoints2, stereoParams1);
[wpsettrial,newPointIndices] = addWorldPoints(wpsettrial,worldPTS);
wpsettrial = addCorrespondences(wpsettrial,1,newPointIndices,indexPairs(:,1));
wpsettrial = addCorrespondences(wpsettrial,2,newPointIndices,indexPairs(:,2));

for h =1: sz(1)

for w = 1:sz(2)
point2track = [pointsFront(h,w,1), pointsFront(h,w,2); pointsRearLeft(h,w,1), pointsRearLeft(h,w,2)];

IDs = [1, 2];
track = pointTrack(IDs,point2track);
triang3D = triangulate([pointsFront(h,w,1), pointsFront(h,w,2)], [pointsRearLeft(h,w,1), pointsRearLeft(h,w,2)], stereoParams1);
% wpSet = addWorldPoints(wpSet,triang3D);
% wpSet = addCorrespondences(wpSet,1,m,m);
% wpSet = addCorrespondences(wpSet,2,m,m);
pointIn3D = [pointIn3D;triang3D];
m = m+1;

end

end
pointIn3D = reshape3D(pointIn3D, board_shape);
%xyzPoints = reshape3D(pointIn3D,board_shape);
s.frontANDleft.PT3D.(frame_i) = pointIn3D;
%s.frontANDleft.PT3DSBA.(frame_i) = xyzPoints;
end
[wpSetRefined,vSetRefined,pointIndex] = bundleAdjustment(wpsettrial,camPoseVSet,[1,3,2],intrinsicsArray, FixedViewIDs=[1,3,2], …
Solver="preconditioned-conjugate-gradient")
end

function [img_name, ptsUS,pts, worldpoints] = reformatData(img_name, pts, board_shape, worldpoints)
%method taken from acinoset code
img_name = img_name(1:strfind(img_name,’ ‘)-1);
img_name = replace(img_name, ‘.’,’_’);
ptsUS = pts;
pts = pagetranspose(reshape(pts, [board_shape, 2]));
pts = pagetranspose(reshape(pts, [board_shape, 2])); %repetition is purposeful

worldpoints = pagetranspose(reshape(worldpoints, [board_shape,2]));
worldpoints = pagetranspose(reshape(worldpoints, [board_shape,2]));
end

function pts = reshape3D(points3D, board_shape)
pts = pagetranspose(reshape(points3D, [board_shape, 3]));
pts = pagetranspose(reshape(pts, [board_shape, 3])); %repetition is purposeful

end Hi everyone! I am working on a project where I need to calibrate multiple cameras observing a scene, to ultimately be able to get 3D points of an object in later videos collected by the same cameras. The cameras are stationary. Importantly, I need to be able to triangulate the checkerboard points from the calibration and then do sparse bundle adjustment on these points to improve the acuracy of the camera pose estimation and 3D checkerboard points from the calibration. Sparse bundle adjustment (bundleAdjustment) can take in either pointTrack objects or worldpointset objects.
I have two calibration sessions (front camera and rear right, and the front camera and rear left – they are in a triangular config) from which I load the stereoParams, I have also stored the useful data in a structure called ‘s’.
I then get the 3D coordinates of the checkerboards, and using the ‘worldpointset’ and feature matching approach. I have included all code (including the code I used to save important variables).
The error I get with the bundleAdjustment function is the following:
Error using vision.internal.bundleAdjust.validateAndParseInputs
The number of feature points in view 1 must be greater than or equal to 51.
Error in vision.internal.bundleAdjust.sparseBA (line 39)
vision.internal.bundleAdjust.validateAndParseInputs(optimType, mfilename, varargin{:});
Error in bundleAdjustment (line 10)
vision.internal.bundleAdjust.sparseBA(‘full’, mfilename, varargin{:});

When I investigated using pointTrack, it seems that it is best used for tracking a point through multiple frames in a video, but not great for my application where I want to track one point through 3 different views at once.

AT LAST–> MY QUESTION:
Am I using worldpointset correctly for this application, and if so, can someone please help me figure out where this error in the feature points is coming from?
If not, would pointTrack be better for my application if I change the dimensionality of my problem? If pointTrack is better, I would need to track a point through the frames of each camera and somehow correlate and triangulate points that way.
**Note, with structure ‘s’ containing many images it was too large to upload (even when compressed), so I uploaded a screenshot of the structure. But hopefully my code helps with context. The visualisation runs though!
load("params.mat","params")
intrinsics1 = params.cam1.Intrinsics;
intrinsics2 = params.cam2.Intrinsics;
intrinsics3 = params.cam3.Intrinsics;
intrinsics4 = params.cam4.Intrinsics;

intrinsicsFront = intrinsics2;
intrinsicsRLeft = intrinsics3;
intrinsicsRRight = intrinsics4;

%% Visualise cameras
load("stereoParams1.mat")
load("stereoParams2.mat")

figure; showExtrinsics(stereoParams1, ‘CameraCentric’)
hold on;
showExtrinsics(stereoParams2, ‘CameraCentric’);
hold off;

%initialise camera 1 pose as at 0, with no rotation
front_absolute_pose = rigidtform3d([0 0 0], [0 0 0]);

%% Get 3D Points
load("s_struct.mat","s")

board_shape = s.front.board_shape;
camPoseVSet = imageviewset;
camPoseVSet = addView(camPoseVSet,1,front_absolute_pose);
camPoseVSet = addView(camPoseVSet,2,stereoParams1.PoseCamera2);
camPoseVSet = addView(camPoseVSet,3,stereoParams2.PoseCamera2);
camposes = poses(camPoseVSet);
intrinsicsArray = [intrinsicsFront, intrinsicsRRight, intrinsicsRLeft];
frames = fieldnames(s.front.points);
framesRearRight = fieldnames(s.rearright.points);
framesRearLeft = fieldnames(s.rearleft.points);
wpSet =worldpointset;
wpsettrial = worldpointset;
for i =1:length(frames) %for frames in front
frame_i = frames{i};
pointsFront = s.front.points.(frame_i);
pointsFrontUS = s.front.unshapedPoints.(frame_i);
container = contains(framesRearRight, frame_i);

j = 1;
if ~isempty(container) && any(container)
pointsRearRight = s.rearright.points.(frame_i);
pointsRearRightUS = s.rearright.unshapedPoints.(frame_i);
pointIn3D = [];
pointIn3Dnew = [];
[features1, validPts1] = extractFeatures(im2gray(s.front.imageFile.(frame_i)), pointsFrontUS);
[features2, validPts2] = extractFeatures(im2gray(s.rearright.imageFile.(frame_i)), pointsRearRightUS);
indexPairs = matchFeatures(features1,features2);
matchedPoints1 = validPts1(indexPairs(:,1),:);
matchedPoints2 = validPts2(indexPairs(:,2),:);
worldPTS = triangulate(matchedPoints1, matchedPoints2, stereoParams2);
[wpsettrial,newPointIndices] = addWorldPoints(wpsettrial,worldPTS);
wpsettrial = addCorrespondences(wpsettrial,1,newPointIndices,indexPairs(:,1));
wpsettrial = addCorrespondences(wpsettrial,3,newPointIndices,indexPairs(:,2));
sz = size(s.front.points.(frame_i));

for h =1: sz(1)

for w = 1:sz(2)
point2track = [pointsFront(h,w,1), pointsFront(h,w,2); pointsRearRight(h,w,1), pointsRearRight(h,w,2)];
IDs = [1, 3];
track = pointTrack(IDs,point2track);
triang3D = triangulate([pointsFront(h,w,1), pointsFront(h,w,2)], [pointsRearRight(h,w,1), pointsRearRight(h,w,2)], stereoParams1);
% [wpSet,newPointIndices] = addWorldPoints(wpSet,triang3D);
% wpSet = addCorrespondences(wpSet,1,j,j);
% wpSet = addCorrespondences(wpSet,3,j,j);
pointIn3D = [pointIn3D;triang3D];
j=j+1;

end

end

pointIn3D = reshape3D(pointIn3D, board_shape);
%xyzPoints = reshape3D(pointIn3D,board_shape);
s.frontANDright.PT3D.(frame_i) = pointIn3D;
%s.frontANDright.PT3DSBA.(frame_i) = xyzPoints;
end
container = contains(framesRearLeft, frame_i);
m=1;
if ~isempty(container) && any(container)

pointsRearLeft = s.rearleft.points.(frame_i);
pointsRearLeftUS = s.rearleft.unshapedPoints.(frame_i);
pointIn3D = [];
pointIn3Dnew = [];
sz = size(s.front.points.(frame_i));
[features1, validPts1] = extractFeatures(im2gray(s.front.imageFile.(frame_i)), pointsFrontUS);
[features2, validPts2] = extractFeatures(im2gray(s.rearleft.imageFile.(frame_i)), pointsRearLeftUS);
indexPairs = matchFeatures(features1,features2);
matchedPoints1 = validPts1(indexPairs(:,1),:);
matchedPoints2 = validPts2(indexPairs(:,2),:);
worldPTS = triangulate(matchedPoints1, matchedPoints2, stereoParams1);
[wpsettrial,newPointIndices] = addWorldPoints(wpsettrial,worldPTS);
wpsettrial = addCorrespondences(wpsettrial,1,newPointIndices,indexPairs(:,1));
wpsettrial = addCorrespondences(wpsettrial,2,newPointIndices,indexPairs(:,2));

for h =1: sz(1)

for w = 1:sz(2)
point2track = [pointsFront(h,w,1), pointsFront(h,w,2); pointsRearLeft(h,w,1), pointsRearLeft(h,w,2)];

IDs = [1, 2];
track = pointTrack(IDs,point2track);
triang3D = triangulate([pointsFront(h,w,1), pointsFront(h,w,2)], [pointsRearLeft(h,w,1), pointsRearLeft(h,w,2)], stereoParams1);
% wpSet = addWorldPoints(wpSet,triang3D);
% wpSet = addCorrespondences(wpSet,1,m,m);
% wpSet = addCorrespondences(wpSet,2,m,m);
pointIn3D = [pointIn3D;triang3D];
m = m+1;

end

end
pointIn3D = reshape3D(pointIn3D, board_shape);
%xyzPoints = reshape3D(pointIn3D,board_shape);
s.frontANDleft.PT3D.(frame_i) = pointIn3D;
%s.frontANDleft.PT3DSBA.(frame_i) = xyzPoints;
end
[wpSetRefined,vSetRefined,pointIndex] = bundleAdjustment(wpsettrial,camPoseVSet,[1,3,2],intrinsicsArray, FixedViewIDs=[1,3,2], …
Solver="preconditioned-conjugate-gradient")
end

function [img_name, ptsUS,pts, worldpoints] = reformatData(img_name, pts, board_shape, worldpoints)
%method taken from acinoset code
img_name = img_name(1:strfind(img_name,’ ‘)-1);
img_name = replace(img_name, ‘.’,’_’);
ptsUS = pts;
pts = pagetranspose(reshape(pts, [board_shape, 2]));
pts = pagetranspose(reshape(pts, [board_shape, 2])); %repetition is purposeful

worldpoints = pagetranspose(reshape(worldpoints, [board_shape,2]));
worldpoints = pagetranspose(reshape(worldpoints, [board_shape,2]));
end

function pts = reshape3D(points3D, board_shape)
pts = pagetranspose(reshape(points3D, [board_shape, 3]));
pts = pagetranspose(reshape(pts, [board_shape, 3])); %repetition is purposeful

end image processing, computer vision, calibration, 3d, worldpointset, bundleadjustment MATLAB Answers — New Questions

​

Trying to translate a piece of Python code below to Matlab. The Python code is provided by Enrichr (see maayanlab.cloud/Enrichr/help#api) as an example for calling its REST APIs.
%{
% – this is the Python code
import json
import requests

ENRICHR_URL = ‘https://maayanlab.cloud/Enrichr/addList’
genes_str = ‘n’.join([
‘PHF14’, ‘RBM3’, ‘MSL1’, ‘PHF21A’, ‘ARL10’, ‘INSR’, ‘JADE2’, ‘P2RX7’,
‘LINC00662’, ‘CCDC101’, ‘PPM1B’, ‘KANSL1L’, ‘CRYZL1’, ‘ANAPC16’, ‘TMCC1’,
‘CDH8’, ‘RBM11’, ‘CNPY2’, ‘HSPA1L’, ‘CUL2’, ‘PLBD2’, ‘LARP7’, ‘TECPR2’,
‘ZNF302’, ‘CUX1’, ‘MOB2’, ‘CYTH2’, ‘SEC22C’, ‘EIF4E3’, ‘ROBO2’,
‘ADAMTS9-AS2’, ‘CXXC1’, ‘LINC01314’, ‘ATF7’, ‘ATP5F1’
])
description = ‘Example gene list’
payload = {
‘list’: (None, genes_str),
‘description’: (None, description)
}
response = requests.post(ENRICHR_URL, files=payload)
%}

% – this is the Matlab code
ENRICHR_URL = ‘https://maayanlab.cloud/Enrichr/addList’;

genes = {‘PHF14’, ‘RBM3’, ‘MSL1’, ‘PHF21A’, ‘ARL10’, ‘INSR’, ‘JADE2’, ‘P2RX7’, …
‘LINC00662’, ‘CCDC101’, ‘PPM1B’, ‘KANSL1L’, ‘CRYZL1’, ‘ANAPC16’, ‘TMCC1’, …
‘CDH8’, ‘RBM11’, ‘CNPY2’, ‘HSPA1L’, ‘CUL2’, ‘PLBD2’, ‘LARP7’, ‘TECPR2’, …
‘ZNF302’, ‘CUX1’, ‘MOB2’, ‘CYTH2’, ‘SEC22C’, ‘EIF4E3’, ‘ROBO2’, …
‘ADAMTS9-AS2’, ‘CXXC1’, ‘LINC01314’, ‘ATF7’, ‘ATP5F1’};

genes_str = strjoin(genes, newline);

description = ‘Example gene list’;

options = weboptions(‘MediaType’, ‘application/json’);

payload(1).(‘list’) = [];
payload(2).(‘list’) = genes_str;
payload(1).(‘description’) = [];
payload(2).(‘description’) = description;

%payload = struct(‘list’, {string(missing), genes_str}, …
% ‘description’, {string(missing), description});

response = webwrite(ENRICHR_URL, payload, options)
But, the translated Matlab code does not work. Any suggestions would be greatly apprecaited.Trying to translate a piece of Python code below to Matlab. The Python code is provided by Enrichr (see maayanlab.cloud/Enrichr/help#api) as an example for calling its REST APIs.
%{
% – this is the Python code
import json
import requests

ENRICHR_URL = ‘https://maayanlab.cloud/Enrichr/addList’
genes_str = ‘n’.join([
‘PHF14’, ‘RBM3’, ‘MSL1’, ‘PHF21A’, ‘ARL10’, ‘INSR’, ‘JADE2’, ‘P2RX7’,
‘LINC00662’, ‘CCDC101’, ‘PPM1B’, ‘KANSL1L’, ‘CRYZL1’, ‘ANAPC16’, ‘TMCC1’,
‘CDH8’, ‘RBM11’, ‘CNPY2’, ‘HSPA1L’, ‘CUL2’, ‘PLBD2’, ‘LARP7’, ‘TECPR2’,
‘ZNF302’, ‘CUX1’, ‘MOB2’, ‘CYTH2’, ‘SEC22C’, ‘EIF4E3’, ‘ROBO2’,
‘ADAMTS9-AS2’, ‘CXXC1’, ‘LINC01314’, ‘ATF7’, ‘ATP5F1’
])
description = ‘Example gene list’
payload = {
‘list’: (None, genes_str),
‘description’: (None, description)
}
response = requests.post(ENRICHR_URL, files=payload)
%}

% – this is the Matlab code
ENRICHR_URL = ‘https://maayanlab.cloud/Enrichr/addList’;

genes = {‘PHF14’, ‘RBM3’, ‘MSL1’, ‘PHF21A’, ‘ARL10’, ‘INSR’, ‘JADE2’, ‘P2RX7’, …
‘LINC00662’, ‘CCDC101’, ‘PPM1B’, ‘KANSL1L’, ‘CRYZL1’, ‘ANAPC16’, ‘TMCC1’, …
‘CDH8’, ‘RBM11’, ‘CNPY2’, ‘HSPA1L’, ‘CUL2’, ‘PLBD2’, ‘LARP7’, ‘TECPR2’, …
‘ZNF302’, ‘CUX1’, ‘MOB2’, ‘CYTH2’, ‘SEC22C’, ‘EIF4E3’, ‘ROBO2’, …
‘ADAMTS9-AS2’, ‘CXXC1’, ‘LINC01314’, ‘ATF7’, ‘ATP5F1’};

genes_str = strjoin(genes, newline);

description = ‘Example gene list’;

options = weboptions(‘MediaType’, ‘application/json’);

payload(1).(‘list’) = [];
payload(2).(‘list’) = genes_str;
payload(1).(‘description’) = [];
payload(2).(‘description’) = description;

%payload = struct(‘list’, {string(missing), genes_str}, …
% ‘description’, {string(missing), description});

response = webwrite(ENRICHR_URL, payload, options)
But, the translated Matlab code does not work. Any suggestions would be greatly apprecaited. Trying to translate a piece of Python code below to Matlab. The Python code is provided by Enrichr (see maayanlab.cloud/Enrichr/help#api) as an example for calling its REST APIs.
%{
% – this is the Python code
import json
import requests

ENRICHR_URL = ‘https://maayanlab.cloud/Enrichr/addList’
genes_str = ‘n’.join([
‘PHF14’, ‘RBM3’, ‘MSL1’, ‘PHF21A’, ‘ARL10’, ‘INSR’, ‘JADE2’, ‘P2RX7’,
‘LINC00662’, ‘CCDC101’, ‘PPM1B’, ‘KANSL1L’, ‘CRYZL1’, ‘ANAPC16’, ‘TMCC1’,
‘CDH8’, ‘RBM11’, ‘CNPY2’, ‘HSPA1L’, ‘CUL2’, ‘PLBD2’, ‘LARP7’, ‘TECPR2’,
‘ZNF302’, ‘CUX1’, ‘MOB2’, ‘CYTH2’, ‘SEC22C’, ‘EIF4E3’, ‘ROBO2’,
‘ADAMTS9-AS2’, ‘CXXC1’, ‘LINC01314’, ‘ATF7’, ‘ATP5F1’
])
description = ‘Example gene list’
payload = {
‘list’: (None, genes_str),
‘description’: (None, description)
}
response = requests.post(ENRICHR_URL, files=payload)
%}

% – this is the Matlab code
ENRICHR_URL = ‘https://maayanlab.cloud/Enrichr/addList’;

genes = {‘PHF14’, ‘RBM3’, ‘MSL1’, ‘PHF21A’, ‘ARL10’, ‘INSR’, ‘JADE2’, ‘P2RX7’, …
‘LINC00662’, ‘CCDC101’, ‘PPM1B’, ‘KANSL1L’, ‘CRYZL1’, ‘ANAPC16’, ‘TMCC1’, …
‘CDH8’, ‘RBM11’, ‘CNPY2’, ‘HSPA1L’, ‘CUL2’, ‘PLBD2’, ‘LARP7’, ‘TECPR2’, …
‘ZNF302’, ‘CUX1’, ‘MOB2’, ‘CYTH2’, ‘SEC22C’, ‘EIF4E3’, ‘ROBO2’, …
‘ADAMTS9-AS2’, ‘CXXC1’, ‘LINC01314’, ‘ATF7’, ‘ATP5F1’};

genes_str = strjoin(genes, newline);

description = ‘Example gene list’;

options = weboptions(‘MediaType’, ‘application/json’);

payload(1).(‘list’) = [];
payload(2).(‘list’) = genes_str;
payload(1).(‘description’) = [];
payload(2).(‘description’) = description;

%payload = struct(‘list’, {string(missing), genes_str}, …
% ‘description’, {string(missing), description});

response = webwrite(ENRICHR_URL, payload, options)
But, the translated Matlab code does not work. Any suggestions would be greatly apprecaited. webwrite rest api, enrichr MATLAB Answers — New Questions

​

I am trying to do my temp_dat variable, have my data in there and my other four variables but doesn’t seem to let me want to do the temp_dat variable. Thank you for your help
Unable to read the entire file. You may need to specify a different format, delimiter, or number of header lines.
Note: readtable detected the following parameters:
‘Delimiter’, ‘,’, ‘HeaderLines’, 0I am trying to do my temp_dat variable, have my data in there and my other four variables but doesn’t seem to let me want to do the temp_dat variable. Thank you for your help
Unable to read the entire file. You may need to specify a different format, delimiter, or number of header lines.
Note: readtable detected the following parameters:
‘Delimiter’, ‘,’, ‘HeaderLines’, 0 I am trying to do my temp_dat variable, have my data in there and my other four variables but doesn’t seem to let me want to do the temp_dat variable. Thank you for your help
Unable to read the entire file. You may need to specify a different format, delimiter, or number of header lines.
Note: readtable detected the following parameters:
‘Delimiter’, ‘,’, ‘HeaderLines’, 0 temp_dat MATLAB Answers — New Questions

​