2014 dxdy logo

Научный форум dxdy

Математика, Физика, Computer Science, Machine Learning, LaTeX, Механика и Техника, Химия,
Биология и Медицина, Экономика и Финансовая Математика, Гуманитарные науки




Начать новую тему Ответить на тему
 
 как переписать код в Matlabe 7
Сообщение13.04.2006, 07:55 
Аватара пользователя


13/04/06
4
Kazan
Ребята у меня проблема с использованием функции fmincon. Раньше в матлабе 6.0 я использовала функцию constr, чтобы минимизировать функционал, зависящий от x0 при ограничениях в виде равенства и неравенства, которые тоже зависят от x0. Ограничения у меня имеют интегральный вид. Подскажите, пожалуйста, могу прислать код.

 Профиль  
                  
 
 
Сообщение13.04.2006, 10:09 
Экс-модератор
Аватара пользователя


23/12/05
12046
А код сильно большой? Может Вы его тут выставите, используя тег code. Кстати, в версии 6.5 присутствуют обе функции, но в описании constr() сказано, что она замещается функцией fmincon() и в последующих версиях уйдет, поэтому предпочтительнее использовать fmincon(). Поскольку у меня стоит 6.5, то я смогу попробовать поэкспериментировать с обеими функциями - давайте, может что и получится.

 Профиль  
                  
 
 
Сообщение18.04.2006, 09:40 
Аватара пользователя


13/04/06
4
Kazan
пожалуйста посмотрите, буду очень признательна, если Вы мне поможите.
%****main**********
clear all;
close all;

global n g m lb1 be cc c2 vmax r0 ep dg ZE ZK consta constb constf0 chetchik
pi2=pi*2;
ka=1.4; % dlya vozduha
hg2=(ka+1)/(ka-1);
n=601;
dg=pi2/(n-1);
g=0:dg:pi2;
ep=2;
%beta=10;
be=0.1; %beta*pi/180;
vmax=1;
vinf=1.0;
% lb=0.6
% mb=sqrt(2/(ka+1))*lb/(sqrt(1-(lb^2/hg2)))
mb=0.5;
r0=0.0;
i=sqrt(-1);
pip=pi/2;
ZE=exp(i*g);
ZK=r0*exp(i*(pi+2*be));
%Kochin-Loicans turbulent PS
constf0=-2.00;
consta=1.17;
constb=4.75;

if mb>0;
c2=0.296
lb=sqrt((ka+1)/2)*mb/sqrt(1+(ka-1)*mb*mb/2) %%???
lb1=2*lb/(1+sqrt(1+4*c2*lb*lb));
else
lb=1.0
lb1=1.0
c2=0.0
end
cc=c2*lb1^2/(1+c2*lb1^2);
options(1)=0;
options(2)=1e-4;
options(3)=1e-4;
%options(14)=1e+10;

m=16;
for k=1:m;
x0(k*2-1)=0.0;
x0(k*2)=0.0;
end

x1=fmincon('fun',x0,options);
%************

function [f,f1]=fun(x0)

global n g m lb1 be cc c2 vmax r0 ep dg ZE ZK consta constb constf0 chetchik %constg0 constb %

i=sqrt(-1);
pip=pi/2;

%generate HI,p,q

HI=log(lb1)+((1-ep)-4*i*exp(i*be)*sin(be)*cc-ZK)./ZE-log(1-ZK./ZE);

for k=1:m;
r=(k+1)^(0);
HI=HI+r*(x0(2*k-1)+i*x0(2*k))./(ZE.^(k+1));
end

p=real(HI);
q=imag(HI);

% calculate
ds=exp(-p).*(2*sin((g)/2)).^(ep-1);
ds=ds-c2*exp(p).*(2*sin(g/2)).^(3-ep).*(2*cos(g/2-be)).^2;

dx=ds.*cos(-q+pip+g-(g-pi)*(ep-1)/2);
dy=ds.*sin(-q+pip+g-(g-pi)*(ep-1)/2);

x(1)=0.0;
y(1)=0.0;
s(1)=0.0;

for k=2:length(g),
s(k)=s(k-1)+(ds(k)+ds(k-1)).*(dg/2);
x(k)=x(k-1)+(dx(k)+dx(k-1)).*(dg/2);
y(k)=y(k-1)+(dy(k)+dy(k-1)).*(dg/2);
end;
v=exp(p).*(2*cos(g/2-be)).*(2*sin(g/2)).^(2-ep);
v=v./(1-c2*v.*v);
vm=max(abs(v));

z=x+i*y;
per=max(s);
%**********************
[cmin,imin]=min(abs(v(3:n-1)));
ista=imin;
vv=abs(v);
ssta=s(ista);

form(ista)=0;
for k=ista-1:-1:1, % на верхней стороне
form(k)=form(k+1)+(vv(k)^(constb-1)+vv(k+1)^(constb-1))*abs((s(k)-s(k+1)))/2;
end;

form1(ista)=0;
for k=ista+1:n, % на нижней поверхности
form1(k)=form1(k-1)+abs((s(k)-s(k-1)))*(vv(k)^(constb-1)+vv(k-1)^(constb-1))/2;
end;
% derivativ on upper side ***************
dv(ista)=(v(ista-1)-v(ista))/abs((s(ista-1)-s(ista)));
for k=ista-1:-1:2;
dv(k)=(v(k-1)-v(k+1))/abs(s(k-1)-s(k+1));
end
dv(1)=(v(1)-v(2))/abs(s(1)-s(2));

% derivativ on bottom side*****************
dv1(ista)=(v(ista+1)-v(ista))/abs(s(ista+1)-s(ista));
for k=ista+1:n-1;
dv1(k)=(v(k+1)-v(k-1))/abs(s(k+1)-s(k-1));
end
dv1(n)=(v(n)-v(n-1))/abs(s(n)-s(n-1));

% %****************************************

% form-param. on uper side

form(ista)=consta/constb;
for k=ista-1:-1:1
form(k)=consta*dv(k)*form(k)/(vv(k)^constb);
end
% %****************************************

% form-parametr on bottom side

form1(ista)=consta/constb;
for k=ista+1:n,
form1(k)=consta*dv1(k)*abs(form1(k))/(vv(k)^constb);
end
fp=form1;

%**********************************************

%**********************


f=max(s);
f1(1)=abs(z(n));
f1(2)=abs(vm-vmax);
f1(3)=constf0-min(form(ista:-1:2));
f1(4)=constf0-min(fp)

 Профиль  
                  
 
 
Сообщение18.04.2006, 10:31 
Аватара пользователя


13/04/06
4
Kazan
.......... ПОГОДИ, КОД НЕМНОГО НЕ ВЕРНЫЙ, НЕ ТО ВЫСТАВИЛА, чуть позже вышлю, спасибо

 Профиль  
                  
 
 
Сообщение19.04.2006, 10:55 
Экс-модератор
Аватара пользователя


23/12/05
12046
А если Вас устраивала функция constr(), почему Вы не можете просто перетащить ее из старой версии - она не built-in, так что это легко должно быть.

 Профиль  
                  
 
 
Сообщение19.04.2006, 11:03 
Экс-модератор
Аватара пользователя


23/12/05
12046
А еще в help-e к Matlab 6.5 нашел образец перехода от функции constr() к fmincon(). Приведу этот раздел help-a полностью

Цитата:
Example of Converting from constr to fmincon
Old Call to constr
Код:
OPTIONS = foptions;

OPTIONS(13) = 2;             % Two equality constraints
OPTIONS(1) = 1;
OPTIONS(9) = 1;
A1 = [ 1 4 -3]; b1 = 2;
A2 = [ 2 5 0]; b2 = 9;
x0 = [1; .5; .8];
LB = []; UB = [];
[X,OPTIONS,LAMBDA,HESS] = ...
   constr('myfuncon',x0,OPTIONS,LB,UB,'mygradcon',A1,b1,A2,b2);

% myfuncon.m
[F, C] = myfuncon(x,A1,b1,A2,b2)
F = x(1) + 0.0009*x(2)^3 + sin(x(3));

C(1,1) = A1*x-b;             % Equality linear constraint
C(2,1) = 3*x(1)^2-1;         % Equality nonlinear constraint

C(3,1) = A2*x-b2;            % Inequality linear constraint
C(4,1) = 7*sin(x(2))-1;      % Inequality nonlinear constraint

% mygradcon.m
[G, DC] = mygradcon(x,alpha)
G = [1;                       % Gradient of the objective
     3*0.0009*x(2)^2;
     cos(x(3))];

DC(:,1) = A1';                % Gradient of the constraints
DC(:,2) = [6*x(1); 0; 0];
DC(:,3) = A2';
DC(:,4) = [0; 7*cos(x(2)); 0];

Цитата:
New Call to fmincon

Код:
OPTIONS = optimset(...
   'Display', 'iter', ...
   'GradCheck', 'on', ... % Check gradients.
   'GradObj', 'on', ...   % Gradient of objective is provided.
   'GradConstr', 'on');   % Gradient of constraints is provided.

A1 = [ 1 4 -3]; b1 = 2;   % Linear equalities
A2 = [ 2 5 0]; b2 = 9;    % Linear inequalities

x0 = [1; .5; .8];
LB = []; UB = [];
[X,FVAL,EXITFLAG,OUTPUT,LAMBDA,GRAD,HESSIAN] = ...
   fmincon(@myfun,x0,A2,b2,A1,b1,LB,UB,@mycon,OPTIONS);

% myfun.m
function [F,G] = myfun(x)
F = x(1) + 0.0009*x(2)^3 + sin(x(3));
G = [1;
     3*0.0009*x(2)^2;
     cos(x(3))];

% mycon.m
function [C,Ceq,DC,DCeq]= mycon(x)
Ceq(1,1) = 3*x(1)^2-1;         % Equality nonlinear constraint
C(1,1) = 7*sin(x(2))-1;        % Inequality nonlinear constraint
DCeq(:,1) = [6*x(1); 0; 0];    % Gradient of equality
                               %    nonlinear constraint
DC(:,1) = [0; 7*cos(x(2)); 0]; % Gradient of inequality
                               %    nonlinear constraint

 Профиль  
                  
 
 
Сообщение19.04.2006, 11:37 
Аватара пользователя


13/04/06
4
Kazan
спасибо, но это я уже изучила. Оказывается проблема в не написании, а в том, что ограничения сложные слишком... СПАСИБО PHOTON ТЕБЕ ЗА ВНИМАНИЕ.

 Профиль  
                  
Показать сообщения за:  Поле сортировки  
Начать новую тему Ответить на тему  [ Сообщений: 7 ] 

Модераторы: Karan, Toucan, PAV, maxal, Супермодераторы



Кто сейчас на конференции

Сейчас этот форум просматривают: нет зарегистрированных пользователей


Вы не можете начинать темы
Вы не можете отвечать на сообщения
Вы не можете редактировать свои сообщения
Вы не можете удалять свои сообщения
Вы не можете добавлять вложения

Найти:
Powered by phpBB © 2000, 2002, 2005, 2007 phpBB Group