function with myblackbox using fminunc

2019 8 30
1 个回答
更新时间：2019 8 30
6 次查看（30 天）

0 个投票

Hello guys! I got a Function F(y(x)) = sum (( yref-y(x))^2) and x(1) = q and x(2)=r and x=[q;r] and yref=0. I wanted to code this function to be used in a multi-objective optimization etc.
my initial idea is:
function F = myblackbox(x)
q = x(1)
r = x(2);
yref = 0;
y = solvemyoptimizationproblem(q,r);
F = somefunctionofy(y);
but i don't know how to use fminunc here to do a blackbox optimization and how to replace those things to have F(y(x)) correctly.

请先登录，再回答此问题。

回答（1 个）

Matt J
Matt J 2019-8-30
编辑：Matt J 2019-8-30

0 个投票

lsqnonlin would be better suited to this,
x0=[q_guess,r_guess];
x=lsqnonlin( @(x) yfunction(x(1),x(2))-yref, x0);

11 个评论
显示 9更早的评论 隐藏 9更早的评论

Ali Esmaeilpour
Ali Esmaeilpour 2019-8-30
would you please write complete function or code? because when I want to run your answer it says undefined 'yofx'
Matt J
Matt J 2019-8-30
The code is as complete as your description allows. In your problem statement
F(y(x)) = sum (( yref-y(x))^2)
You haven't told us what y(x) looks like.
Ali Esmaeilpour
Ali Esmaeilpour 2019-8-30
y is closed-loop response and x=[q;r] and q and r are weighted matrices.
Ali Esmaeilpour
Ali Esmaeilpour 2019-8-30
my closed-loop response contains a 1 by 100 matrix of numerical values and x is a decision variable here and x=[q;r] and I'm going to tune those q and r
Matt J
Matt J 2019-8-30
编辑：Matt J 2019-8-30
Then yfunction, in my post, should be the function that takes x as input and returns the corresponding 1x100 response.
Ali Esmaeilpour
Ali Esmaeilpour 2019-8-30
how can I find a function for my closed-loop response y? I got closed-loop response:
y = [-2,-2.06844254709663,-2.10111830954110,-2.10207250524360,-2.07526182090699,-2.02451336009587,-1.95348936029213,-1.86566104239566,-1.76428609501489,-1.65239283602609,-1.53276900719044,-1.40795510853059,-1.28024204134041,-1.15167219215486,-1.02404365938112,-0.898917175767882,-0.777625266347119,-0.661283291728917,-0.550801952158784,-0.446900918848503,-0.350123272116316,-0.260850455112571,-0.179317482440462,-0.105628171524862,-0.0397701935839232,0.0183702310931665,0.0689941424455793,0.112376278061898,0.148852238329798,0.178806457148514,0.202661035909770,0.220865478976649,0.233887353198917,0.242203878667164,0.246294444722945,0.246634033866838,0.243687526483402,0.237904851041483,0.229716937775089,0.219532428513973,0.207735091324267,0.194681885774748,0.180701622852288,0.166094162756435,0.151130093813326,0.136050836531366,0.121069118234489,0.106369765666474,0.0921107653476128,0.0784245442470077,0.0654194263502996,0.0531812239571451,0.0417749259245079,0.0312464485645638,0.0216244182313353,0.0129219582599530,0.00513845621642096,-0.00173870934149472,-0.00773250023679152,-0.0128746228204809,-0.0172040168372426,-0.0207654382076332,-0.0236081426698606,-0.0257846749601092,-0.0273497662907451,-0.0283593408257317,-0.0288696309906560,-0.0289363993725932,-0.0286142641583750,-0.0279561239979665,-0.0270126774854485,-0.0258320316788283,-0.0244593937199454,-0.0229368391498630,-0.0213031503020680,-0.0195937183993203,-0.0178405024207615,-0.0160720382595013,-0.0143134917503068,-0.0125867493914005,-0.0109105408257964,-0.00930058753949017,-0.00776977247844266,-0.00632832582476125,-0.00498402235945812,-0.00374238646270093,-0.00260690103919371,-0.00157921715665090,-0.000659361544098163,0.000154060515616681,0.000863668978079722,0.00147312127494728,0.00198693439542896,0.00241031792821887,0.00274901888833613,0.00300917888601873,0.00319720397995884,0.00331964734521793,0.00338310467267475,0.00339412212095973];
Ali Esmaeilpour
Ali Esmaeilpour 2019-8-30
I give you the code that I got this close-loop response from and mux1 adn mux2 are closed-loop responses I mean y=mux1 in the following program:
clc;
clear;
close all;
%% Time structure
t0 = 0;
tf = 10;
dt = 0.1;
time = t0:dt:tf-dt;
%% System structure ===> x(k+1) = A*x(k) + B*u(k) + G*w(k)
A = [1.02 -0.1;0.1 0.98];
B = [0.5 0;0.05 0.5];
G = [0.3 0;0 0.3];
[nx,nu] = size(B);
%% MPC parameters
Q = eye(nx);
R = eye(nu)*50;
Qf = eye(nu)*50;
N = 26;
%% Problem parameters
M = blkdiag(Q,R);
w1 = randn(1,100);
w1 = w1-mean(w1);
muw1 = mean(w1);
sigmaw1 = var(w1);
w2 = randn(1,100);
w2 = w2-mean(w2);
muw2 = mean(w2);
sigmaw2 = var(w2);
w = [w1;w2];
sigmaw = [sigmaw1 0 ;0 sigmaw2];
x = randn(2,1);
mux{1} = [-2;2];
sigmax{1} = [1,0;0,1];
h1x = [-1/sqrt(5);-2/sqrt(5)];
h2u = [-0.4/sqrt(1.16);1/sqrt(1.16)];
epsilon = 0.5;
g1 = 3;
g2 = 0.2;
c = 0;
%% Main loop
for t = t0:dt:tf-dt
c = c+1;
sigmax = sdpvar(repmat(nx,1,N+1),repmat(nu,1,N+1));
p = sdpvar(repmat(4,1,N-1),repmat(4,1,N-1));
pN = sdpvar(2,2);
F = sdpvar(repmat(2,1,N+1),repmat(2,1,N+1));
K = sdpvar(repmat(2,1,N+1),repmat(2,1,N+1));
muu = sdpvar(repmat(nu,1,N),repmat(1,1,N));
sigmau = sdpvar(repmat(2,1,N),repmat(2,1,N));
constraint = [];
objective = 0;
for k = 1:N-1
mux{k+1} = A*mux{k} + B*muu{k};
objective = objective + 1/2*(trace(M*p{k}));
constraint2 = [sigmax{k+1} (A*sigmax{k}+B*F{k}) G*sigmaw;(A*sigmax{k}+B*F{k})' sigmax{k} zeros(2);(G*sigmaw)' zeros(2) sigmaw]>=0;
constraint3 = [sigmau{k} (F{k});(F{k})' sigmax{k}]>=0;
constraint4 = h1x'*mux{k}<=(1-(0.5*epsilon))*g1-((0.95/2*epsilon*g1*(1-0.95))*h1x'*sigmax{k}*h1x);
constraint5 = h2u'*muu{k}<=(1-(0.5*epsilon))*g2-((0.95/2*epsilon*g2*(1-0.95))*h2u'*sigmau{k}*h2u);
constraint6 = [p{k} [sigmax{k};F{k}] [mux{k};muu{k}];[sigmax{k};F{k}]' sigmax{k} zeros(2,1);[mux{k};muu{k}]' zeros(1,2) ones(1)]>=0;
constraint7 = [F{k} K{k};sigmax{k} eye(2)];
constraint = [constraint,constraint2,constraint3,constraint4,constraint5,constraint6,constraint7];
end
objective = objective + (1/2*trace((trace(Q*pN))));
constraint8 = h1x'*mux{N}<=(1-(0.5*epsilon))*g1-((0.95/2*epsilon*g1*(1-0.95))*h1x'*sigmax{N}*h1x);
constraint9 = [pN-sigmax{N} mux{N};mux{N}' 1];
constraint = [constraint,constraint8,constraint9];
option = sdpsettings('solver','sedumi');
a = optimize (constraint,objective,option);
u = value(muu{1});
U1(c) = u(1,1);
U2(c) = u(2,1);
U3 = [U1;U2];
mux3 = value(mux{1});
mux1(c) = mux3(1,1);
mux2(c) = mux3(2,1);
F = value(F{1});
X = value(sigmax{1});
f1{c} = F;
X1{c} = X;
K = F*inv(X);
% K=value(K{1})
K1{c} = K;
K2 = value(K1{c});
x = value(x);
u1 = u+K2*(x-mux3);
u2{c} = u1;
x = A*x+B*u1+G*w(:,c);
mux{1} = value(mux{2});
sigmax{1} = value(sigmax{2});
ud = value(u2{c});
ud1(c) = ud(1,1);
ud2(c) = ud(2,1);
x1(c) = x(1,1);
x2(c) = x(2,1);
end
%% Plot results
figure(1)
subplot(4,1,1)
plot(time,ud1,'LineWidth',2)
xlabel('time')
ylabel('ud1')
title('Input signal')
grid on
subplot(4,1,2)
plot(time,ud2,'LineWidth',2)
xlabel('time')
ylabel('ud2')
grid on
subplot(4,1,3)
plot(time,x1,'LineWidth',2)
xlabel('time')
ylabel('x1')
title('State 1')
grid on
subplot(4,1,4)
plot(time,x2,'LineWidth',2)
xlabel('time')
ylabel('x2')
title('State 2')
grid on
figure(2)
subplot(4,1,1)
plot(time,x1,'r',time,mux1,'g','LineWidth',2)
xlabel('time')
legend({'x1','Mean of x1'},'Location','northeast')
grid on
subplot(4,1,2)
plot(time,x2,'r',time,mux2,'g','LineWidth',2)
xlabel('time')
legend({'x2','Mean of x2'},'Location','northeast')
grid on
subplot(4,1,3)
plot(mux1,mux2,'LineWidth',2)
title('State space phase plane')
xlabel('mux1')
ylabel('mux2')
grid on
subplot(4,1,4)
plot(U1,U2,'LineWidth',2)
title('Control space phase plane')
xlabel('U1')
ylabel('U2')
grid on
figure
subplot(4,1,1)
plot(time,mux1,'b','LineWidth',2)
title('Mean of first state')
xlabel('Time')
ylabel('mux1')
grid on
subplot(4,1,2)
plot(time,mux2,'b','LineWidth',2)
title('Mean of second state')
xlabel('Time')
ylabel('mux2')
grid on
subplot(4,1,3)
plot(time,U1,'r','LineWidth',2)
title('First input')
xlabel('Time')
ylabel('U1')
grid on
subplot(4,1,4)
plot(time,U2,'r','LineWidth',2)
title('Second input')
xlabel('Time')
ylabel('U2')
grid on
Ali Esmaeilpour
Ali Esmaeilpour 2019-8-30
for this code you should have YALMIP in your matlab path and also sedumi as solver installed
Matt J
Matt J 2019-8-30
编辑：Matt J 2019-8-30
how can I find a function for my closed-loop response y?
To use lsqnonlin or any other Matlab optimizer, all you need is code that generates y for a given choice of [x(1), x(2)], which you seem to have already. For lsqnonlin, you also need this function to be differentiable, and it is not clear from the complexity of your code if that is the case.
SInce we are not sure it is differentiable, then you could use fminsearch instead, which doesn't use derivatives, and works well for a small number of unknowns,
x=fminsearch(@(x) norm(yfunction(x)-yref), [q_guess, r_guess])
Ali Esmaeilpour
Ali Esmaeilpour 2019-8-30
so I put that fminsearch at the end of my main code?
Matt J
Matt J 2019-8-30
fminsearch will search for the optimal x. You put it wherever you need the optimization to occur.

请先登录，再进行评论。

类别

在 帮助中心 和 File Exchange 中查找有关 Problem-Based Optimization Setup 的更多信息

产品

版本

R2018b

标签

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by